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A new numerical framework for solving conservation aws is be- 
ing developed. This new framework differs substantially m both 
concept and methodology from the well-established methods ue 
finite difference, finite volume, finite element, and spectral methods^ 

It is conceptually simple and designed to overcome several y 
limitations of the above traditional methods. A two-level scheme 
for solving the convection-diffusion equation 

r'JU/rif + a - n i) 2 u/i>X 2 = 0 {/X > 0) 

is constructed and used to illuminate major differences between the 
present method and those mentioned above. This exp//c/f scheme, 
referred to as the a-p scheme, has two independent marching van 
abies u” and (u,)? which are the numerical analogues of u and 
Mai U, n). respectively. The a-p scheme has the unusual prop- 
erty that its stability is limited only by the CFL condition, i.e., it 
independent of p. Also it can be shown that the amplification factors 
of the a-p scheme are identical to those of the Leapfrog scheme ,f 
u = 0 and to those of the DuFort-Frankel scheme if a - 0. These 
coincidences are unexpected because the a-p scheme and the above 
classical schemes are derived from completely different perspec- 
tives, and the a-p scheme does not reduce to the above class ‘” 
schemes in the limiting cases. The a-p scheme is extended to solve 

the ID time-dependent Navier-Stokes equations of a perfect gaS; 

Stability of this explicit solver also is limited only by e 
tion. In spite of the fact that it does not use (.) any techn '^® 
related to the high-resolution upwind methods, and (it) any ad ho 
parameter, the current Navier-Stokes solver is capaMe i o 9 a "® ra 

ing highly accurate shock tube solutions. Particularly, for high-Rey 

olds-number flows, shock discontinuities can be resolved with n 
one mesh interval. The inviscid {p = 0) a-p scheme is reversib n 
time. It also is neutrally stable, i.e., free from numerical dissipat on 
Such a scheme generally cannot be extended to soNe the Euler 
equations. Thus, the inviscid version is modified.. Stability of this 
modified scheme, referred to as the scheme, is limited by the 
CFL condition andOstSl, where a is a special parameter that 
controls numerical dissipation. Moreover, if « - 0, the amplification 
factors of the a-, scheme are identical to those of the Leapfrog 
scheme, which has no numerical dissipation. On the other hand, if 
K = 1 the two amplification factors of the a-e scheme become 
the same function of the Courant number and the phase angle. 
Unexpectedly, this function also is the amplification fact0 '' 0 
highly diffusive Lax scheme. Note that, because the Lax scheme 
very diffusive and it uses a mesh that is staggered in time, a two- 
level scheme using such a mesh is often associated with a highly 


ffusive scheme. The a-c scheme, which also uses a mesh stag 
sring in time, demonstrates that it can also be a scheme with no 
jmerical dissipation. The Euler extension of the a-c scheme 
ability conditions similar to those of the a-e. scheme itself. It has 
,e unusual property that numerical dissipation at all mesh poi 
an be controlled by a set of local parameters. Moreover, it is capable 
f generating accurate shock tube solutions with the CFL number 
anainq from Close to 1 to 0.022 « 1995 Academic Press, Inc. 


1. INTRODUCTION 

The method of space-time conservation element and solution 
element [1-11] is a new numerical framework tor solving 
conservation laws. This new approach differs substantially in 
both concept and methodology from the well-established meth- 
ods. i.e.. finite difference, finite volume, finite element, and 
spectral methods [ 12-16). It is conceived and designed to over- 
come several key limitations of the above traditional methods 
Thus we shall begin this paper with a discussion of severa 
considerations that motivate the current development: 

(a) A set of physical conservation laws is a collection of 
statements of flux conservation in space-time. Mathematically, 
these laws are represented by a set of integral equations. The 
differential form of these laws is obtained from the Integra 
form with the assumption that the physical solution is smooth. 
For a physical solution in a region of rapid change (e.g.. a 
boundary layer), this smoothness assumption is difficult to real- 
ize by a numerical approximation that can use only a limite 
number of discrete variables. This difficulty becomes even 
worse in the presence of discontinuities (e.g., shocks). Thus, a 
method designed to obtain numerical solutions to the differen- 
tial form without enforcing flux conservation is at a fundamental 
disadvantage in modeling physical phenomena with high-gradi- 
ent regions. Particularly, it may not be used to solve flow 
problems involving shocks. Contrarily. a numerical solution 
obtained from a method that also enforces flux-conservation 
locally (i.e., down to a computational cell) and globally (i.e., 
over the entire computational domain) will always retain the 
basic physical reality of flux conservation even in a region 
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involving discontinuities. For this reason, the enforcement of 
both heal amt global flux conservation in space and time is a 
tenet in the current development. To meet this requirement, 
first we define a set o f solution elements which are subdomains 
in the space-time computation domain. Within each solution 
element, any physical flux vector is then approximated in terms 
°f Nome s ' m P le smooth functions. In the last step, we divide 
the computational domain into conservation elements and de- 
mand that any flux be conserved over any space-time region 
that is the union of any combination of these elements. Note 
that a solution element generally is not a conservation element 
and vice versa. 

Among the traditional methods, finite difference, finite ele- 
ment. and spectral methods are designed to solve the differential 
form of the conservation laws. Note that the set of integral 
equations usually solved in a finite-element scheme is equiva- 
lent to the differential form of the conservation laws assuming 
certain smoothness conditions. However, these integral equa- 
tions generally are different from the integral equations repre- 
senting the conservation laws. Even if they are cast into a 
conservative form, the resulting flux-conservation conditions 
generally do not represent the physical conservation laws. 

The Unite volume method is the only traditional method 
designed to enforce flux conservation. A finite-volume scheme 
may enforce flux conservation in space only, or in both space 
and time. As a preliminary' to this enforcement, a flux must be 
assigned at any interface separating two neighboring conserva- 
tion cells. In a typical finite-volume scheme, it is evaluated by 
extrapolating or interpolating the mesh values at the neigh- 
boring cells. This evaluation generally requires an ad hoc choice 
of a special flux model among many models available ( 17-19], 
Generally numerical results obtained are dependent on which 
model one chooses. Also this process of interpolation and ex- 
trapolation generally is time consuming and has some undesir- 
able side effects which will be discussed shortly. 

Contraiily, by defining conservation elements wisely and 
considering the the spatial derivatives of dynamic variables as 
independent variables, current flux evaluation at an interface 
is carried out without interpolation or extrapolation. It is an 
integral part of the solution procedure. 

(b) Space and time traditionally are treated separately in 
the time marching schemes. Generally one obtains a system of 
ordinary differential equations with time being the independent 
variable after a spatial discretization. As an example, elements 
in the finite element method usually are used for spatial discreti- 
zation. These elements are domains in space only. 

Because flux conservation is fundamentally a property in 
space-time, space and time are unified and treated on the same 
looting in the present method. Thus, conservation elements and 
solution elements used in the time-dependent version of the 
present method are domains in space-time. The significance of 
this unified approach cannot be overemphasized. As will be 
shown, it makes it easier for a numerical analogue to share the 
same space-time symmetry of the physical laws. 


{ c> ln a finite-difference scheme, derivatives at mesh points 
are expressed in terms of mesh values of dependent variables by 
using finite-difference approximations. The accuracy of these 
approximations, especially those of higher-order accuracy, gen- 
erally is excellent as long as dependent variables vary slowly 
across a mesh interval. However, it may not be adequate if 
these variables vary too rapidly. Thus, in a high-gradient region, 
e.g„ a boundary layer, accuracy may demand the use of an 
extremely fine mesh. In turn, a prohibitively high computing 
cost may result. 

The present method avoids the above pitfall by (i) expressing 
the numerical solution within a solution element as an expansion 
in terms of certain base functions, and (ii) considering the 
expansion coefficients as the independent numerical variables 
to be solved for simultaneously. For simplicity, Taylor’s expan- 
sions will be used in the present paper. For this special case, 
the expansion coefficients are interpreted as the numerical ana- 
logues of the derivatives. Note that (i) van Leer (20| also has 
attempted to improve accuracy by introducing two independent 
numerical variables for each independent physical variable, and 
GO the current solution procedure has no resemblance with 
those used in compact difference schemes. 

(d) The numerical variables used in a spectral method, i.e., 
the expansion coefficients, are global parameters pertaining to 
the entire computational domain. As a result, a spectral method 
generally (i) lacks local flexibility and thus may be applied 
only to problems with simple geometry, and (ii) is hindered by 
the fact that it must deal with a full matrix that is difficult 
to invert. 

By design, only local parameters will be used in the present 
method. Moreover, solution elements and conservation ele- 
ments are defined such that the set of discrete variables in 
any one of the numerical equations to be solved generally is 
associated with only two neighboring solution elements. The 
exception to this general rule occurs only in the situation in 
which numerical dissipation is introduced deliberately. Even 
in this special case, only the discrete variables associated with 
a few immediately neighboring solution elements will enter any 
equation to be solved. Thus, a scheme developed using the 
present method generally has the simplest stencil and one needs 
only to deal with a very sparse matrix if the scheme is implicit. 
Moreover, the maximum number of solution elements involved 
in a numerical equation of the current discretization framework 
is independent of the order of accuracy of a particular scheme. 
The order of accuracy can be raised by using a Taylor’s expan- 
sion of higher order as the approximated solution within a 
solution element. Contrarily, the order of accuracy of a classical 
finite-difference scheme generally can be increased only by 
using variables of more mesh points in each of its equations. 
Usually, a side effect of this practice is an increase in numerical 
dissipation, a subject to be discussed shortly. Also it may be 
difficult to implement a high-order finite-difference scheme 
near a boundary because there are no real mesh points outside 
the boundary. The above discussions also point to another im- 



CONSERVATION-SOLUTION ELEMENT 


297 


portant advantage of the present method, i.e., the specification 
of initial/houndary conditions generally is simpler , more flexi- 
ble, and more accurate than that associated with a traditional 
method . It is simpler because a smaller stencil is used 1 6). 
Furthermore, it is more flexible and accurate because the spatial 
derivatives of dynamical variables, which are considered as 
independent numerical variables in the present method, can 
now be specified directly. Note that the current emphasis in 
reducing the size of the stencil is also consistent with a funda- 
mental physical reality, i.e., in the absence of body force, direct 
physical interaction occurs only among the immediate 
neighbors 

(e) With a few exceptions, numerical dissipation generally 
appears in a numerical solution of a time-marching problem. 

In other words, the numerical solution dissipates faster than the 
corresponding physical solution. For a nearly inviscid problem, 
e.g., flow with a high Reynolds number, this could be very 
serious because numerical dissipation may overwhelm physical 
dissipation and cause a complete distortion of solutions. One 
may argue that numerical dissipation can be reduced by increas- 
ing the order of accuracy of the scheme used. However, because 
the order of accuracy of a scheme is generally determined with 
the aid of Taylor’s expansion, and the latter is valid only for 
a smooth solution, it has meaning only for a smooth solution. 
Thus the use of a scheme of higher-order accuracy may not 
reduce numerical dissipation associated with high-frequency 
Fourier components of a numerical solution. This is the reason 
that the Leapfrog scheme, w hich is tree from numerical dissipa- 
tion, can outperform schemes with higher-order accuiacy in 
solving some wave equations [21]. 

In a study of finite- difference analogues of a simple convec- 
tion equation [2|, it was shown that a numerical analogue will 
be free from numerical dissipation if it does not violate certain 
space-time invariant properties of the convection equation. In 
other w r ords, numerical dissipation may be considered as a 
result of symmetrv-breaking by the numerical scheme. Because 
of its intrinsic nature of space-time unity, the current framework 
is an excellent vehicle for constructing a numerical analogue 
that shares the same space-time invariant properties with the 
physical equation. 

It is recognized that a certain amount of numerical dissipation 
may be needed to prevent large dispersive errors [22] that are 
often caused by the presence of high-frequency disturbances 
(such as round-off errors). Therefore, in the present paper we 
shall construct a model scheme for a simple convection equation 
in which its numerical dissipation is controlled by a single 
adjustable parameter. The numerical dissipation is shut off 
when this parameter is set to zero. Furthermore, an Euler solver 
will be constructed such that its numerical dissipation at all 
mesh points can be controlled by a set of local pararnetet s. 

(f ) High-resolution upwind methods 116] form a special 
class of the finite volume method. In these methods, the flux 
at an interface separating two neighboring conservation cells 


is also evaluated using a process of interpolation and extrapola- 
tion. This process generally is heavily dependent on characteris- 
tics-based techniques. For the ID time-dependent case, the 
characteristics are curves in space-time, and the coefficient 
matrix associated with the Euler equations |23] also can be 
diagonalized easily. As a result, these techniques are easy to 
apply. However, for multidimensional cases, the characteiistics 
are 2D or 3D surfaces in space-time [24|. Moreover, the coeffi- 
cient matrices cannot be diagonalized simultaneously by the 
same matrix [23], Because of the above complexites, applica- 
tion of these techniques to multidimensional problems is much 
more diffucult. Furthermore, high-resolution methods generally 
require the use of ad hoc parameters, e.g., flux-limiters and / 
or slope-limiters, and other ad hoc techniques. These ad hoc 
techniques may lead to numerical dissipation which varies trom 
one place to another and from one Fourier component to an- 
other. In other words, numerical solutions may suffer annihila- 
tion of sharply different degrees at different locations and differ- 
ent frequencies [5, 25]. Also, these techniques generally are 
also difficult to apply in a space of higher dimension. 

Although only the ID time-marching schemes are con- 
structed in the present paper, the current framework is devel- 
oped to solve multidimensional problems. In order that ID 
schemes can be extended to become multidimensional schemes 
in a straightforward manner, simplicity and generality weigh 
heavily in the development of the present method. Thus, we 
do not use characteristics-based techniques, and also try to avoid 
using ad hoc techniques. Note that, except the Navier- .Stokes 
solver, other ID schemes described in the present paper have 
been extended to become their 2D counterparts [7, 8) (the 
extension of the Navier-Stokes solver will be dealt with in a 
separate paper). Also, because of the similarity in their design, 
each of the 2D schemes described in [7, 8] shares with its ID 
version virtually the same fundamental characteristics. Further- 
more, it is shown in [71 that a 2D Euler time-marching solver, 
which uses a uniform stationary mesh, is capable of generating 
highly accurate solutions lor a 2D shock reflection problem 
used by Helen Yee and others [26f Specifically, both the inci- 
dent and the reflected shocks can be resolved by a single 
data point without the presence of numerical oscillations near 
the discontinuity. 

In addition to being difficult to apply in a space of higher 
dimension, the concept of characteristics generally is also not 
applicable to the Navier-Stokes equations, which is non-hyper- 
bolic in nature. Therefore, the decision not to use characteris- 
tics-based techniques also makes it easier for the present frame- 
work to solve the Navier-Stokes equations. 

This completes the discussion of the motivation for the cur- 
rent development. In summary, the development is guided by 
the following requirements: (i) to enforce both local and global 
flux consevation in space and time with flux evaluation at an 
interface being an integral part of the solution procedure and 
requiring no interpolation or extrapolation; (ii) space and time 
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are unified and treated on the same footing; (iii) mesh values 
of dependent variables and their derivatives are considered as 
independent variables to be solved for simultaneously; (iv) to 
use only local discrete variables; (v) solution elements and 
conservation elements be defined such that the simplest stencil 
will result; (vi) to minimize numerical dissipation, a numerical 
analogue should be constructed, as much as possible, to be 
compatible with the space-time invariant properties of the corre- 
sponding physical equations; and (vii) to exclude the use of 
the characteristics-based techniques, and to avoid the use of ad 
hoc techniques as much as possible. It is the purpose of this 
paper and its follow-ups [6—8] to show’ that the above require- 
ments can be met with a simple unified numerical framework. 

For any reader who is interested in getting an advance idea 
on how ; simple the present method can be, he is referred to the 
computer program listed at the end of the present paper. It 
is a shock-tube-problem solver constructed using the present 
method. The simplicity of the solver is easily appreciated by 
a comparison of the listed program and a typical program 
associated with high-resolution upw ind methods Not only is 
the listed program much smaller in size (it is self-contained 
and the main loop contains only 33 lines), hut it contains no 
Fortran statements such as * 'if ' ’ ‘ ‘amax, ' ’ and ' ‘am in ' ’ which 
are used so often in the programs implementing high-resolution 
methods . The absence of the above Fortran statements in the 
listed program results from the effort in avoiding the use of 
the ad hoc techniques in the development of the present method. 
In spite of its simplicity, it will be shown in Section 7 that the 
present solver is capable of generating highly accurate shock 
tube solutions. 


2. THE a-p SCHEME 

In this section, we consider a dimensionless form of the ID 
convection— diffusion equation, i.e.. 


(ht (hi d~u 

— + u - r . ft — = 0, 

dr ax ^ hx- 


( 2 . 1 ) 


where the convection velocity a . and the viscosity coefficient 
M (— 0) are constants. Let X] = a , and x 2 = / be considered as 
the coordinates of a two-dimensional Euclidean space E 2 . By 
using Gauss' divergence theorem in the space-time E 2 , it can 
be shown that Eq. (2.1) is the differential form of the integral 
conservation law 


^ ^ hits = 0. (2.2) 

As depicted in Fig. 1, here (i) S(V) is the boundary of an 
arbitrary space-time region V 7 in £\, (ii) h = (an - fjihu/hx, //) 
is a current density vector in E : , and (iii) ds = da n with da 
and n, respectively, being the area and the outward unit normal 
of a surface element on S(V). Note that (i) h * ds is the space- 


?=(x, t) 



FIG. I. A surface element ds and a line segment dr on the boundary 5(V) 
of a volume V in a spacc-timc E : . 


tune flux of h leaving the region V through the surface element 
ds, and (ii) all mathematical operations can be carried out 
as though E : were an ordinary two-dimensional Euclidean 
space. 

At this juncture, note that the conservation law' given in Eq. 
(2.2) is formulated in a form in which space and time are unified 
and treated on the same footing. This unity of space and time 
is also a tenet in the following numerical development. It is a 
key characteristic that distinguishes the present method from 
most of the traditional methods . 

Let fl denote the set of mesh points (j\ n) in E 2 (dots in Fig. 
2(a)), where n = 0, ±i, ±1, ±f, ±2, ±£ and, for each w, 

j = n — L n ± f, n ± f There is a solution element (SE) 

associated with each (;, n) G II. Let the solution element SE (f 
n) be the interior of the space-time region bounded by a dashed 
curve depicted in Fig. 2(b). It includes a horizontal line segment, 
a vertical line segment, and their immediate neighborhood. For 
the following discussions, the exact size of this neighborhood 
does not matter. 

For any (x, t) G SE(y, n), u(x. /), and h(.r, t), respectively, 
are approximated by m*(.v, t\ j, n) and h*(.v, rj, n) which we 
shall define shortly. Let 

u*{x, r, j, n) = u’j + (u,)'}(x - xj) + ( u,)"j(t - i"), (2.3) 

where (i) u”, f«J", and («,); are constants in SE (j, n ). and (ii) 

( Xj , /") are the coordinates of the mesh point (j, n). Note that 


u*(. Xj, t"; j , n) = u i. 


„ ftw*(.v, r, j, n) 


<)x 

'>u*(x, k j, n) 
at 


= (u,)% 


= («,)?. 


(2.4) 


Moreover, if we identify «;,(«,)", and («,);. respectively, with 
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j-i i J+i 



n + 1 
n + 1/2 
n 

n -1/2 
n - 1 



(j,n) 


(J-1/2, 
n - 1/2) 


(j,n) 

A. D 


a .□ 


(i + 1/2, 

n - 1/2) 


of the diffusion term. Combining Eqs. (2.3) and (2.5), 
one has 


u*(x, t: j, n) = + (M,)"|(jr - .v f ) - a(t - /“)), 

(,v. /) G SE(./', ii). 


( 2 . 6 ) 


Because h = (cut — ftdu/hx, u), we dehne 

h*(x, V, j. n) = (au*(x y /; j, n) ^ 

- ilHu*(x, f, j, n)/cix, u* (x, /; j, it)). 

Let E : be divided into nonoverlapping rectangular regions 
(see Fig. 2(a)) referred to as conservation elements (CEs). As 
depicted in Figs. 2(c) and 2(d), the CE with its top-right (top- 
left) vertex being the mesh point (j, it) G 11 is denoted by 
CE it) (CE,(j, //)). Obviously the boundary of CE (j, n) 
(CE. ( /. n)), excluding two isolated points B and C (C and D). 
is formed by the subsets of SE( /. n) and SE( J — n — i) 
(SE( j + i it - i)). The current approximation of Eq. (2.2) is 


e 


(j- 1 / 2 , 

n + 1/2) 



(i,n) 


(j + 1/2, 
n + 1/2) 


FIG. 2. The SEs and CEs of type I: (a) The relative positions' of SEs and 
CEs; (b) SE( j. «>: (c) CE (/, n): (d) CE,(./. n); (e) CK.(; - i, n + i)'. (f) 
CE ( j + i n + i). 


the values of u , du/()x, and I'luA'lt at (x,. /"), the expression on 
the right side of Eq. (2.3) becomes the first-order Taylor’s 
expansion of u(x, 1) at (x r l' r ). As a result of these considera- 
tions, (it,)", and ( u,)1 will be considered as the numerical 
analogues of the values of u, c)u/c)x. and clu/dt at (.v, , /"), 
respectively. 

We shall require that u = u*(x. r, j , n) satisfy Eq. (2.1) 
within SE(y, it). As a result of Eq. (2.4), this implies that 

(«,)'/ = -«(«,)"• (2 - 5 > 

Because Eq. (2.3) is a first-order Taylor's expansion, the diffu- 
sion term in Eq. (2.1) has no counterpart in Eq. (2.5). As a 
result, the diffusion term has no impact on how w*(.v, f, j. n) 
varies with time within SE(j, it). However, as will be shown 
shortly, through its role in the numerical analogue of Eq. 
(2.2), it does influence time-dependence of numerical solutions. 
Note that, for a higher-order scheme, how n*(.v, /; j, it) varies 
with time within SE (j, it) will be influenced by the presence 


F\(/ii)=T h* • ds = 0 (2.8) 

J ,Y((T , (.;,«) I 

for all ( /, n) £ H. In other words, the total flux leaving the 
boundary of any conservation element is zero. Note that the 
flux at any interface separating two neighboring CEs is calcu- 
lated using the information from a single SE. As an example, 
the interface AC depicted in Figs. 2(c) and 2(d) is a subset of 
SE (j, n). Thus the flux at this interface is calculated using the 
information associated with SE(/, n). Also note that an SE is 
the interior of a space-time region. Thus the vertices B , C, and 
D, strictly speaking, do not belong to any SE. As a result, h* 
is not defined at these points. However, contributions to the 
above integral from these isolated points are zero no matter 
what values of h* are assigned to them. For this reason, one 
may simply exclude them from the above surface integration. 

Because the surface integration across any interface separat- 
ing two neighboring CEs is evaluated using the information 
from a single SE, obviously the local conservation condition 
Eq. (2.8) will lead to a global conservation relation, i.e„ the 
total flux leaving the boundary of any space-tune region that 
is the union of any combination of CEs will also vanish . 

Because each S(CEAj, n)) is a simple closed curve in E 2 
(see Fig. 1 ), the surface integration in Eq. (2.8) can be converted 
into a line integration. Let 

g* = (-*/*, an* - /jL(lu*/bx ), dr = (dx, dt). (2.9) 

Thus, dr is normal to ds and points in the tangential direction 
of the line segment joining the two points (xJ) and (.v + dx , 
t + dt). Because ds = ±(du -dx) | L p. 14], w^e have 
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h* els = ±g* -t/r, (2. 10) 

where the upper (lower) sign should be chosen if the 90° rotation 
from ils to elr is in the counterclockwise (clockwise) direction. 
By combining Eqs. (2.8) and (2.10), one concludes that 

= g*-dr. (2.11) 

J .SIC L . l/.ni) 


Note that the notation c.c. indicates that the line integration 
should be carried out in the counterclockwise direction. Substi- 
tuting Eq. (2.6) into Eq. (2.11), and using the fact that the 
boundary of a CE is formed by the subsets of two SEs, one has 


(A.v) 


iFAj, n) 


± I T ) K 1 ~ + (1 - ir - $)(u x )?:{j?\ (2.12) 


2(1 n- v) 

+ — — (uj ~ up x jj) % 


\ - V -(1 - V 1 - f) \ 

1 - v 1 -(1 4 v){\ - v 2 ~ £) . (2.17) 

1 - v 1 4 £ I - v 2 4 f / 

Because numerical variables at a higher time level can be 
evaluated in terms of those at a lower time level by using Eq. 

(2.14) , it defines a marching scheme. Furthermore, because 
this scheme models Eq. (2.1) which is characterized by two 
parameters a and /it, hereafter it will be referred to as the 
a-fx scheme. 

As a preliminary for future developments, we apply Eq. 

(2.14) successively and obtain 

q(j,n 4 1 ) = (Q+) 2 q(j ~ 1 ,/?) 

+ (Q+G- + Q (? )q(y>0 (2.18) 

+ (Q ) 2 q(y + 1 ,/i) ( 1 — v 2 + £#<)). 

A result of Eq. (2.18) is 



q<4 w + 1 ) — > q( j, n) asA/— »(), (2.19) 


where 


if cu \x, and A.v are held constant. The proof follows from the 
fact that 


det 
V — 


a At 
Ax 


4pAt_ 

(Ax) 2 ' 


(2.13) 


Note that (i) the parameter v is the Courant number, and (ii) a 
more efficient method of flux evaluation will be presented later 
in this section. 

With the aid of Eqs. (2.8) and (2.12), u) and («,)" can be 
solved in terms of up and (wj;;, 1 /^ if 1 - \r + £ ^ 0; i.e., 
for all SE(y, /f). 


q O '. ») = 0»q( / ~ in - 

i) 


+ Q-*lU + i n 

-b (l - V 1 + £*0). 

(2.14) 

Here 

, del 

( "" ) 


q(y. n) = 


(2.15) 

for all (y, n) G 12, and 

/ 1 4 v 

1 

1 

'JYf 


~ del / 1 \ / 

Q' = U I -d - ") 

I 

1 

1 

T 

(2.16) 

^ \l-r4f 

1 - V 2 + f / 



0, (Q+Q- + 0 00-^1, 
(Q ) 2 ->0 as At — ► 0, 


( 2 . 20 ) 


if r/, /it, and Ax are held constant. 

Alternatively, Eq. (2.19) can be proved using the fact that 
the total flux of h* leaving the boundary of any space-time 
region that is the union of any combination of CEs vanishes. 
Consider the union of CE+(y, n 4 1) and CE.(y + i n 4 £) 
(see Fig. 2). This union is a rectangle with the vertices ( j 4 
h n + 1), (y, n + 1), (y, n) and (y 4 - n). The flux leaving 

this rectangle through its two vertical edges approaches zero 
as At —> 0. Because the total flux leaving its boundary vanishes, 
one concludes that the total flux leaving its two horizontal edges 
also approaches zero as A/— > 0. In other words, the flux entering 
the rectangle through the lower horizontal edge approaches that 
leaving through the upper horizontal edge as At — 0 . Because 
these two fluxes are evaluated using q(y, n) and q (/, /if 1), 
respectively, the above limiting condition implies a limiting 
relation between q(y, n) and q(y, n 4 1 ). Similarly, by consider- 
ing the union of CE (y, n 4 1 ) and CE + (j - J, n 4 £), one 
obtains another limiting relation for q( y. n) and q(y, n 4 1 ). 
Equation (2.19) is a result of the above two limiting relations. 

The a-fx scheme has several nontraditional features. They 
are summarized in the following remarks: 


and 


(a) Space and time are unified and treated on the same footing 
in the construction of the a-fx scheme. 
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(b) The expansion coefficients u] and (u x )] in Eq. (2.6) are 
treated as independent variables; i.e., («,))' is not expressed in 
terms of u'f s by using a finite-difference approximation. 

(c) As a result of Eq. (2. 12), each of the conservation condi- 
tions FAj , n) involves only numerical variables associated with 
two neighboring SEs. This fact remains true for a scheme 
of higher-order accuracy in which Eq. (2.3) is replaced by a 
Taylor s expansion of higher-order. The contrast with the finite 
difference method and its physical significance were dicussed 
in Section 1. 

(d) The a-fx scheme has the simplest stencil, i.e., a triangle 
with a vertex at the upper time level and the other two vertices 
at the lower time level. Equation (2. 14), which relates numerical 
variables at these vertices, was derived using the flux conserva- 
tion conditions F (7, n) = 0. Because the flux at an interface 
separating two neighboring CEs is evaluated using information 
of a single SE, no interpolation or extrapolation is required. 
Moreover, accuracy of flux evaluation is enhanced by requiring 
that a = w*U, /; 7, n) satisfy Eq. (2.1) within SE (/, n). This 
makes the use of character! stics-based techniques less nec- 
essary. 

(e) The a-fx scheme uses a mesh that is staggered in time. 
As will be explained in Appendix A, for a two-level scheme 
using such a mesh, e.g., the Lax scheme [12, p.97], generally 
the numerical variable at (7, n + 1) does not approach that at 
(j n n) as A/ 0. if a, fi, and Ajc are held constant. This is a 
key reason why the Lax scheme is very diffusive when the 
Courant number v is small. According to Eq. (2.19), the a-fx 
scheme is an exception to the above general rule. 

(f) Equation (2.1) can be solved numerically using the 
Leapfrog/DuFort-Frankel scheme [12, p.1611. This scheme is 
reduced to the Leapfrog scheme [12, p.100] if diffusion is 
absent (i.e., fx = 0), and to the D u Fort- Fran kel scheme [12, 
p.l 141 if convection is absent (i.e., a = 0). It is well known 
that a solution of any of the above schemes is tormed by two 
decoupled solutions with each being associated with a mesh 
that is also staggered in time. Traditionally the von Neumann 
stability analysis for the above schemes is performed without 
taking into account this decoupled nature [12]. In Appendix A, 
it is performed separately for each decoupled solution using 
the mesh depicted in Fig. 2(a). It is show ; n that the amplification 
factors of the Leapfrog/DuFort-Frankel scheme are 


marching steps. The reason behind this definition is that the 
mesh points at the time levels n and n + 1 are not staggered. 
Let 1 - v 2 ¥> 0. Then the amplification factors G\ f of the 
current a-fi scheme (see Eq. (6.9)) are identical to those given 
by Eq. (2.21) except that the parameter f should be replaced 
by | = £/( l - p 2 ). Because (i) | = 0 if /x = 0, and 

(ii) v = 0 and thus | = £ if a = 0, one concludes that G'L are 
completely identical to those of the Leapfrog scheme if fx = 

0, and to those of the Du Fort -Fran kel scheme if a = 0. These 
coincidences are unexpected because the a-jx scheme and the 
above classical schemes are derived from completely different 
perspectives. Moreover, the a-jx scheme is a tw ; o-level scheme 
with tw'O variables u'f and (u V )J ! associated wfith the mesh point 
(y, /?), while the above classical schemes are three-level schemes 
with a single variable u” associated with the same point. 

Because the amplification factors of the inviscid a-fx scheme 
(i.e., the a-fx scheme with fx — 0) are identical to those of the 
Leapfrog scheme, the former, as in the case of the latter, is 
neutrally stable (i.e., free of numerical dissipation) if ir < \. 
Note that the case with jx = 0 and v 1 — 1 is ruled out by the 
assumption 1 — v + £ ^ 0 of Eq. (2.14). Similarly, the pure- 
diffusion a-p scheme (i.e., the a-fx scheme with a = 0), as in 
the case of the DuFort-Frankel scheme, is unconditionally 
stable. Furthermore, it is proved in Section 6 that the stability 
of the general ci-fx scheme, as in the case of the Leapfrog/ 
DuFort-Frankel scheme, is independent of fx , and restricted 
only by the CFL condition , i.e., v 2 < 1. The a-p scheme is the 
only two-level explicit scheme know n to the author to possesss 
the above properties. Also it will be shown later that the same 
stability condition is retained by a natural ID time-dependent 
Navier-Stokes extension of the a-p scheme. 

Because stability of the a-p scheme is restricted only by the 
CFL condition, the stability bound for A; is proportional to A*. 
In contrast, the stability condition of a typical classical explicit 
scheme generally is more restrictive than the CFL condition. 
For a small mesh Reynolds number, the stability bound for A t 
is approximately proportional to (A,v)" for the MacCormack 
scheme [12, p. 102]. 

Because a neutrally stable numerical analogue of the pure 
convection equation 

— + a — — 0 ( 2 . 22 ) 

bt bx 


A ± = 



[£cos(0/2) - r>sin(0/2) 


± V[£cos(0/2) - />sin(tf/2)] 2 + 



( 2 . 21 ) 


Here 0, -n < 0 < n [1, p.30], is the phase angle variation per 
A.v. Note that, in the present paper, the amplification factors 
are defined to be those between the time levels n and n + 1, 
i.e., they are the amplification factors of the solution after two 


usually becomes unstable when it is applied to a nonlinear 
inviscid generalization of Eq. (2.22), the inviscid a-fx scheme 
will be modified in Section 3 such that it can be extended 
to model the Euler equations. In this new version, numerical 
dissipation is introduced in a w F ay that allows its magnitude to 
be adjusted by a special parameter. 

(g) The conservation relations for CE_(7 “ b n + 2) and 
CE (7 + i n + 5) (see Figs. 2(e) and 2(f)) are 

Ffj - i /?+£) = 0, FAj + i n + £) = 0, (2.23) 
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respectively. Combining Eqs. (2.12) and (2.23) and assuming 
1 ~ ir ~ £ # 0, one has 




( 1 + v ~ (1 — v 2 + £) 

1 - V 2 "(1 - v)(\ - v 2 + g) 
1 — v 2 — £ 1 — V 2 — £ 


(2.25) 


( 1 - 1 - v 2 + £ \ 

-(1 - V 1 ) -(1 + v)( \ - V 1 + £) . (2.26) 

1 - - f I - v 2 - f / 

Equation (2.24) defines a backward marching scheme, i.e., the 
numerical variables at the time level n are determined in terms 
of those at the time level in + J). Recall that both the forward 
marching scheme Equation (2.14) and the backward marching 
scheme Eq. (2.24) are derived using the same set of conservation 
relations. As a matter of fact, Eqs. (2.14) and (2.24) are equiva- 
lent if (1 - v~y ¥= (£)- is assumed. For the above reason, the 
a-fx scheme may be referred to as a two-way marching scheme. 
For the case fx > 0, it will be proved in Section 6 that the a-/x 
scheme cannot be stable for both the forward and the backward 
marching directions, except for the singular case v 2 = 1 which 
is also on the threshold of instability. Thus, for all practical 
purposes the viscous a-fx scheme is irreversible in time. On 
the other hand, it is neutrally stable for both the forward and 
backward marching directions, and thus is reversible in time, 
\i (x = 0 , and i* 2 < 1. Again, the a-/x scheme is the only two- 
level explicit two-way marching scheme known to the author. 

(h) Several invariant properties of Eq. (2.1) with respect to 
space and time are discussed in [2]. In the same paper, these 
properties are also defined for the numerical analogues of Eq. 
(2. 1 ). It is also shown that the neutral stability of several finite- 
difference analogues of Eq. (2.22) can be established by using 
their invariant properties with respect to space-time inversion. 
Because solutions of Eq. (2.22) do not dissipate with time, it 
is not surprising that solutions of a numerical analogue also 
will not dissipate with time, i.e., the scheme is neutrally stable, 
if it shares with Eq. (2.22) some space-time invariant properties. 

It will be shown in a future paper that the a-fx scheme shares 
with Eq. (2.1) the same space-time invariant properties. Also 
note that these invariant properties are closely linked with the 
other properties discussed in (a), (e), (f), and (g). 

This completes the discussion on nontraditional features of 
the a-fx scheme. In the following, it will be shown that this 
scheme can also be constructed from a completely different 




FIG. 3. The SEs and CEs of type II: (a) 1 he relative positions of SEs and 
CEs: (b) SE( /, n): (c) CEO', «); (d) Three neighboring CEs. 

perspective. As a part of this construction, SEs and CEs of 
different types will be used and discussed. 

In the new construction, the locations of mesh points (dots in 
Fig. 3(a)) are identical to those used in the original construction. 
However. SE (y, n) is defined to be the interior of a rhombus 
centered at (y, n) (see Fig. 3(b)). CE(y, n) is the union of 
SE(y, n) and its boundary. Readers are warned not to confuse 
the sides of the rhombus with the characteristics of Eq. (2.22). 
Any one of these sides is simply a line segment joining two 
points of intersection (not marked by dots) of horizontal and 
vertical mesh lines. For any ( v, t) G SE(y, n ), m(jc, t) and h(x, 
/), respectively, again are approximated by u *(a, ;; j, n) and 
h*(.v, /; j, n) which are defined by Eqs. (2.3) and (2.7), respec- 
tively. However, Eq. (2.5) will be derived from a consideration 
of flux conservation. 

Let Eq. (2.2) be approximated by 

^ h*-tfs = 0, (2.27) 

where V* is the union of any combination of CEs. Because an 
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SE is the interior of a CE, h* is not defined on S(V*), the 
boundary of V*. As a result, the above surface integration is 
to be carried out over a surface that is in the interior of V* 
and immediately adjacent to S(V*). A necessary condition of 
Eq. (2.27) is that, for all (y, n ) £ 12. 

h* ■ ds = 0; (2.28) 

.SVC Hi j f./i n 


the numerical solution uniformly satisfies the differential form 
of the conservation law Eq. (2.2). 

With the aid of Gauss’ divergence theorem, Eq. (2.30) im- 
plies that the surface integration of h* over any closed surface 
located within any SE vanishes. As a result, 

h*</s = (), l h* ■ ds = 0, (2.31) 

S{AAHC) J SiAA'H'C ) 


i.e., the total flux leaving any conservation element is zero. 

Note that the center of a current SE no longer sits on an 
interface separating two CEs. It coincides with the center of a 
CE. Thus h* at one side of an interface is evaluated using 
information from one SE, while that at the other side is evalu- 
ated using information from another SE. As an example, h* at 
BC and B'C depicted in Fig. 3(d), respectively, are evaluated 
using information from SE (y, n) and SE(y —£,/? — \). Another 
necessary condition for Eq. (2.27) is the equality between the 
fluxes entering and leaving any interface. This can be seen by 
applying Eq. (2.27) separately to two neighboring CEs, and then 
to their union. Obviously the local flux conservation relations at 
all interfacs, and within all CEs (i.e., Eq. (2.28)) are equivalent 
to the global conservation relation Eq. (2.27). The equations 
representing the above conservation conditions are the numeri- 
cal equations to be solved. Note that, in the current construction, 
a flux is not preassigned at an interface using an interpolation 
or extrapolation of information from both sides of this interface. 
The present method of interface flux evaluation obviously is 
different from that used in the finite volume method which was 
discussed in Section 1 . 

By using Eqs. (2.3) and (2.7), one concludes that, for any 
(a, t) £ SE(y, «), the divergence of h* in E 2 is 

„ , . jet d[au*(x t /; y, n) — a <lu*(x, /; y, n)Hlx\ 

V * n* — ; 

()X 


r)w*U, t\ y, n) 
in 

= aUi x )J + (u,yj. 


(2.29) 


where the triangles A ABC and A A' B'C' are those depicted in 
Fig. 3(d). Because the net flux of h* entering an interface from 
both sides vanishes, the sum of the flux leaving CE(y, n) through 
BC and that leaving CE (y - - h) through B'C vanishes. 

Thus, Eq. (2.31) implies that F (j. n) = 0, where F.(j , «) is 
defined in Eq. (2.11). Similarly, it can be shown that 
F.(f n ) = 0. 

Assuming Eqs. (2.3) and (2.7), it has been shown that both 
Eqs. (2.5) and (2.8) can be derived using Eq. (2.27). Conversely, 
Eq. (2.27) also follows from Eqs. (2.5) and (2.8). Obviously 
both the forward marching scheme Eq. (2.1 4) and the backward 
marching scheme Eq. (2.22) can also be obtained by assuming 
Eqs. (2.3), (2.7), and (2.27). 

Note that the equivalence between Eq. (2.27) and the pair 
of equations Eqs. (2.5) and (2.8) hinges on the fact that 
V • h* = 0 within an SE of either type I or type II. As will be 
shown immediately, this condition can be used to simplify 
evaluation of the flux across a simple curve that lies entirely 
within an SE of either type. 

According to the top expression given in Eq. (2.29), 
Vh* = 0 implies that there exists a function if/*(x, t\ j , /?) 
such that 


f)t/f(jc, t\ y, n) 

in 


= au*{x, /; y, n) ~ p 


Hu*{x, t\ y, n) 
hx 


(2.32) 


and 


hipjx, t\ j, n) 

in 


m*(jc, t\ y, n) 


(2.33) 


Because (w.,)" and (w f )" are constants within an SE, Eq. (2.29) 
implies that V • h* is also a constant. Thus Eq. (2.28) coupled 
with Gauss’ divergence theorem implies that, within any SE, 

V h* = 0. (2.30) 

Equation (2.5) is a direct result of Eqs. (2.29) and (2.30). 

Note that Eq. (2.30) follows from Eq. (2.28) because «*(*, 
t : ; y, n) dehned in Eq. (2.3) is a hrst-order Taylor’s expansion 
For a higher-order expansion, the condition that Eq. (2.30) 
being valid uniformly within an SE is stronger than Eq. (2.28) 
For the general case, the stronger condition should be imposed. 
Because Eq. (2.30) is the numerical analogue of Eq. (2.1), the 
imposition of the stronger condition ensures that, within an SE, 


for any (x, t) £ SE(y, a?)- Substituting Eq. (2.6) into Eqs. (2.32) 
and (2.33), one concludes that, up to an arbitrary constant, 

(Mj/ r 

<AU, t\ y, n) = y~ {[(.v - Xj ) - a(t - /")!“ 

(2.34) 

+ 2 flit - /")} - u';\(x - Xj) - a(t - m 

Moreover, with the aid of Eq. (2.9), Eqs. (2.32) and (2.33) 
imply that 

g*-</r = #. (2.35) 

Let (jc, t) £ SE(y, n) and (x\ t f ) £ SE(y, n). Let T be a simple 
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FKS. 4. A simple curve T joining U, /) and (.v\ t'). 


curve joining (a\ t) and (.v\ /'), and lying entirely within 
SE (j, n) (see Fig. 4). Then Eqs. (2.10) and (2.35) imply that 

h* • c Is = fi(x\ f\ j. n) — t (f(x % t\ j , «). (2.36) 

Here we assume that ds points to the right of T if one moves 
forward from (a, /) to (a\ /') (see Fig. 4). Equation (2.36) states 
that the flux of h* across the curve F is given by the difference 
in the values of if at its two end-points. For this reason. ijj(x , 
t: j , n) will be referred to as the potential function associated 
with SE (y, n). Obviously, Eq. (2.12) can be obtained using 
Eq. (2.36). 

Note that a generalized a-p scheme with a moving mesh 
was constructed in [ 1 1. This scheme is reduced to the present 
a-jji scheme when the mesh becomes stationary. In [1], the 
generalized scheme is subjected to a thorough theoretical and 
numerical analysis on stability, dissipation, dispersion, consis- 
tency, truncation error, and accuracy. It is shown that it has 
many advantages over the MacCormack and the Leapfrog/ 
Dufort-Frankel schemes. Particularly, by using a new discrete 
Fourier error analysis, it is shown that the generalized scheme 
is more accurate than the Leapfrog/DuFort-Frankel scheme 
by one order (in a sence defined in (1|) in both initial-value 
specification and the main marching scheme. Other key results 
of [ 1 ] are summarized in the following remarks: 

(a) For the generalized scheme, (i) stability and accuracy 
can be improved, and (ii) dissipation and dispersion can be 
reduced, if the space-time mesh is allowed to evolve with the 
physical variables such that the local convective motion of 
physical variables relative to the moving mesh is kept to a 
minimum. 

(b) For a numerical analogue of Eq. (2.22) that has both 
principal and spurious amplification factors, a numerical solu- 
tion with periodic boundary conditions is the sum of a principal 
solution and a spurious solution |1, p.32]. Only the principal 
solution contributes to the accuracy of the scheme. Note that 
(i) the behaviors of the principal and the spurious solutions as 
functions of time are determined by the principal and spurious 
amplification factors, respectively; (ii) both two amplification 
factors of the present inviscid a-fx scheme are of unit magnitude; 
and (iii) given an accurate initial-value specification, the spuri- 


ous solution at / = 0 generally is very small compared with the 
principal solution. As a result, the spurious solution of the 
present a-fi scheme generally is negligible. Furthermore, for 
the inviscid a-jx scheme, it is shown that [ 1, pp. 36-37] (i) the 
principal solution has no dispersion if v — 0 or in the limit of 
1 ; and (ii) each Fourier component of the principal solution 
has a convection velocity not more than a and not less than 
(2ht)a for all phase angles and ir < I. In other words, the 
dispersion associated with the inviscid a-pi scheme is small 
compared with that associated with a typical finite-difference 
scheme. 

In conclusion, a model scheme has been constructed from 
two different perspectives using SEs and CEs of different types. 
Using either perspective, one can say that a numerical solution 
generated using the current framework satisfies (i) the differen- 
tial form of the conservation law uniformly within an SE, and 
(ii) the integral form over any region that is the union of any 
combination of CEs. The second perspective that used the SEs 
and CEs of type II depicted in Fig. 3 was used in the initial 
development of the present method ( 1 ). In addition, it also was 
adopted to develop several new solvers for the 2D steady, 
incompressible Navier-Stokes equations |4, 9-11]. However, 
in these new solvers, the CEs and SEs depicted in Fig. 3 are 
replaced by CEs and SEs of rectangular shape in the 2D spatial 
computational domain. It was shown that, for a laminar channel 
flow with Re ; = 100, an accurate solution can be obtained by 
using as few as six SEs across the channel. 

Because (i) the first perspective is easier to use in constructing 
explicit schemes, and (ii) the schemes to be discussed in the 
present paper are exclusively explicit, the first perspective will 
be adopted in the present paper hereafter. 

3. THE a-e SCHEME 

The inviscid a-p. scheme is neutrally stable and reversible 
in time. It is well known that a neutrally stable numerical 
analogue of Eq. (2.22) generally becomes unstable when it is 
extended to model the Euler equations. It is also obvious that 
a scheme that is reversible in time cannot model a physical 
problem that is irreversible in time, e.g., an inviscid flow prob- 
lem involving shocks. In this section, we assume fx — 0 and 
attempt to modify the inviscid a-p scheme such that it can be 
extended to model the Euler equations. 

The current path of development is almost identical to that 
given in Section 2. We continue to assume Eqs. (2.3)— (2.7), 
and use SEs of type I depicted in Fig. 2. In addition to f x — 0, 
the only other modification is the replacement of the assumption 
F~(j, n) = 0 by 


E . (7, n) 


8(1 - ir)(kx) 2 
4 


(duff. 


(3.1) 


where 8 is a parameter independent of numerical variables, and 
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(duj; = (i)| + (Mv)/-i l /fl - U'r m ~ «" 


(3.2) 


(d) Because ( m t j'/f is a numerical analogue of du/dx at 
( j ± £, /? — I), the simple average 


In other words, we add two terms of the same magnitude but 
with opposite signs, respectively, to the right sides of the origi- 
nal conservation conditions FA.h n) = 0 and F (y, n) = 0 . 
The beauty of this modification will be fully explained later in 
this section. For now it suffices to say that this modification 
injects a higher-order finite-difference error into the inviscid 
a-jji scheme. It breaks the space-time symmetry of the latter. 
In turn, numerical dissipation is introduced as a result of this 
symmetry breaking. Because the magnitude of the terms added 
in this modification is controlled by e, numerical dissipation 
is controlled by z in the modified scheme just as physical 
dissipation is controlled by p in the a-p scheme. Note that, as 
a result of Eq. (3.1) and the assumption p = 0, the modified 
scheme is characterized by two parameters a and z. Thus, 
hereafter it will be referred to as the a-z scheme. Also note 
that, because there is no upwind bias in the a-z scheme, upwind 
bias is not the source of numerical dissipation. Additional 
remarks on Eqs. (3.1) and (3.2) are: 

(a) By definition, F.(y, n) and F ( j , n) represent total (luxes 
leaving CE + (y, n) and CE J, /?), respectively (see Figs. 2 (c) 
and 2(d)). Because FAj, n) # 0 if CE.(y, ft) and CE_(y, 
n) generally are no longer conservation elements in the 
a-z scheme. 

(b) Let CE(y, n) be the union of CE f (y, n) and CE (y, n) 
(see Fig. 5(b)). Note that this definition of CE(/, n) differs from 
that given in Section 2 and depicted in Fig . 3(c). Let 


(i)| (wjj.i’/f + UOlflr 1 

is a numerical analogue of dufftx at ( /. n — i), the midpoint 
of a line segment joining (y + h n - s) and ( j — i, n — i) 
(see Fig. 2(a)). Note that (y, n — i) <£ (1 if (y, n) G (1. Also 
note that 

(u'filr: ~ u'l ,'n )/A.v 

is a central-difference analogue of hulhx at (y, n - %)■ Thus, 
(du x Y; represents the difference of two numerical analogues of 
hii I fix at the same mesh point (y, n — 4). By using Taylor's 
expansion at (y, n - 5 ), it can be shown that (duffi = 
0 [(Ajr) 2 ], if (mJ-T j 1 /? are identified with dw(.v,ii/ 2 , C l/ 2 )/r).v, 
respectively. Hereafter a quantity is denoted by 0 [(Ajt) , 1 if 
there exists a constant C > 0 such that the absolute value of 
this quantity | A.v | ; for all sufficiently small |A.v|. Note 
that we have constructed an expression of 0 [(A.v) 2 | without 
explicitly introducing the factor ( Ajc) 2 . This natural construction 
leads to the simple stability conditions to be given in Eq. (3.14). 
It is possible only because there are two discrete variables 
w" and {u x yj associated with the mesh point (y, n). 

(e) Equation (3.1) could have been written as F.(y, n) = 
±e\du x ) n j with z 1 =^£(1 - r 2 )(A.v) 2 /4, However, this simplified 
expression would lead to much more complicated equations 
later. 


F( /.«) = <£ h*ds. (3.3) 

Because the net flux entering the interface separating CE< 
(/, n) and CE_ (y, n) is zero. F(y, n) is the sum of FAj, n) 
and F (y, n). With the aid of Eq. (3. 1), we have 

F(y, n) = FA j , n) + FA j, n) = 0; (3.4) 


This completes the discussion of Eqs. (3.1) and (3.2). Now, 
let \ - v 2 ^ 0. Then Eqs. (2.12), (3.1), and (3.2) can be used 
to obtain the current counterparts of Eqs. (2.14) and (2.18). 
They are 


q O',/!) = A/ + q(y - in - i) 

+ M q(y + b n - b) ( 1 — v 1 ^ 0) 


(3.6) 


i.e., the total flux leaving CE(y, n) vanishes. As a result, CE 
(y, n) is a conservation element in the a-z scheme. Note that 
Eq. (3.4) leads to a global conservation relation in the form of 
Eq. (2.27), where V* is the union of any combination of these 
new CEs. 

(c) Because £ = 0 if p = 0, Eq. (3.4) coupled with Eq. 
( 2 . 12 ) implies that 


and 

q(y, n + 1 ) = (M . ) 2 q( j ~ L n) 

+ (A FM + M.MAqifin) (3.7) 

+ (M ) : q(y + L n) ( 1 - v 1 7 ^ 0), 


«; = (a + + (i - 


4 - 


Ax(l » — — I («,)./ -In - («,)"- i'/rl- 


respectively. Here 


(3.5) 


M., = (?) 


1 + V I — v 1 
z — 1 2 e — 1 + v 


Thus, u" is independent of z. 


and 


(3.8) 
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Obviously, M , — Q L if e = 0 and £ = 0. Furthermore, the 
limiting condition given in Eq. (2. 19) is still valid if we assume 
that e = e{At) and lim^ 0 s(\t) = 0. However, unlike the a- 
fx scheme, the a-G scheme is not a two-way marching scheme 
if e 0. 

Equation (3.6) represents a pair of equations. The first is Eq. 
(3.5). With the aid of Eqs. (2.5) and (2.13), the second equation 
can be expressed as 

ui x y; = utjim - «;-i,2)/Ajc + (2c - nutoj;. (3.io) 

Here 

= MjVi'/f + (Stl2)(u,)Wd : (3. 1 1 ) 

i.e., Uj’t\t 2 is a first-order Taylor’s approximation of u at (j ± 
5 , n). Thus, the expression on the right side of Eq. (3. 10) is the 
sum of a central -difference approximation of dufdx at (j, n) 
and the extra term (2c - [)(du x YJ. Because (du x )” = 0[(A.*) 2 |, 
the presence of this extra term will not lower the order of 
accuracy of the entire sum as an approximation of dufdx at 
( j , n ). Also note that this extra term vanishes when g = i 
while the term associated with ( du x ))' in Eq. (3.1) vanishes 
when e — 0. 

Next we shall study the influence of c on the stability and 
numerical dissipation of the as scheme. Let G <2> and G <2) be 
the principal and spurious amplification factors of the as 
scheme, respectively. Then, it will be shown in Section 6 that 

G { r>= [A,(e, (3.12) 


by comparing Eqs. (2.21), (3.12), and (3.13) w-ith £ = 0 and 
e = 0. 

Also, we have 

A±(l, v , 0) = cos(0/2) - />sin(0/2). (3.15) 

Thus, G { f = G i2 J when £ = I . Moreover, it is shown in Appen- 
dix A that the coalesced amplification factor is identical to that 
of the Lax scheme . Note that, like the Leapfrog scheme, a 
solution of the Lax scheme is also composed of two decoupled 
solutions with each being associated with a mesh that is stag- 
gered in time. However, because the Lax scheme is a two-level 
scheme, it does not have a spurious amplification factor. 

Thus, at one extreme, i.e., when e = 0, become the 
amplification factors of the Leapfrog scheme, which is free of 
numerical dissipation. At another extreme, i.e., when e — 1 , 
Gf and G i2i coalesce into one and it becomes the amplification 
factor of the Lax scheme, which is notorious for its large diffu- 
sive errors. From the above observations, one may infer the 
conclusion that will be established shortly, i.e., the a-e scheme 
becomes more diffusive as the value of e increases . Note that , 
because the Lax scheme is very r diffusive and uses a mesh that 
is staggered in time, a two- lev el scheme using such a mesh is 
usually associated with a highly diffusive scheme 1 27 1 . The 
a-G scheme demonstrates that it can also be a scheme with no 
diffusive error! 

As a result of Eq. (3.14), the expression under the radical 
sign in Eq. (3.13) is nonnegative. Thus, it can be shown that 

• - |G'i1 = xAe, v, 9) = e{(l - p 2 ) sin 2 (0/2) 4- 2 cos(0/2) 

X [(1 - e)cos(0/2) (3.16) 

+ V(1 - e)[(l - e)cos 2 (0/2) + (1 - tX) sin’(0/2)|]}. 


Because solutions to the physical equation Eq. (2.22) do not 
dissipate with time, a numerical analogue to Eq. (2.22) is said 
to be free of numerical dissipation if its solutions also do not 
^ I ^ dissipate with time, i.e., its amplification factors are of unit 
magnitude. As a result, numerical dissipation of the a-e scheme 
may be measured by 1 — |G ,2) |, i.e., ^.(s. v, 0). Obviously the 
a-e scheme is free of numerical dissipation if e = 0. Also, by 
using Eqs. (3.14) and (3.16), it is shown in Section 6 that, for 
all 6 with — tt < 6 < 77, and all e and v satisfying Eq. (3.14), 
we have 

are necessary and sufficient conditions for the stability of the 
a-e scheme. Thus, Eq. (3.14) will be assumed in the remainder 
of this section. 

It was pointed out in Section 2 that the amplification factors 
of the Leapfrog scheme are identical to those of the inviscid 
a-(x scheme. Because the latter scheme is a special case of the anc * 
a-G scheme with s = 0, G [ : ] become the amplification factors 
of the Leapfrog scheme when c - 0. This fact can be reverified 


0 < 0) 4- 4e( 1 — e) cos‘(0/2) 

< X (e, v, 0 ) < min{ 1 , 4 e} 


( 3 . 17 ) 


0 ^ 0) — e(l “ v 2 ) sin 2 (0/2). 


with 

A^(e, v , 0) = ecos(0/2) — /Vsin(0/2) 

± V(1 - £ )[ ( 1 — e)cos 2 (0/2) + (1 - ir) sin : (0/2)J. 

Also it will be proved that 

0 < £ < 1 and ir < 1 


(3.18) 
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The significance of Eqs. (3.17) and (3. 1 8) is discussed in the 
following remarks: 

(a) Note that (i) the behaviors of the principal and the spurious 
solutions of the a-e scheme are determined by its principal and 
spurious amplification factors, respectively; and 00 because 
() < s(l - e) if 0 < e < 1, Eq. (3.17) implies that *.(e. v , 

0 ) < x ( e< v , 0). Thus, the spurious solution will not dissipate 
more slowly than the principal solution. Let e be not too close 
to 0 or 1. Then Eq. (3.17) also implies that the Fourier compo- 
nents of the spurious solution with smaller | 0 | i.e., longer wave- 
length. will dissipate much faster than those of the principal 
solution. In other words, the spurious solution will rapidly 
disappear from the long- wavelength components of a numerical 
solution. Note that * (i v, 0) = 1. Thus, the long-wavelength 
components of the spurious solution are annihilated almost 
completely in a single time step if e = i i.e., if the last term 
in Eq. (3.10) is dropped. 

(b) The upper bound of * + (b. v , d) given in Eq. (3.18) is 
proportional to sin : (0/2). As a result, the long-wavelength Fou- 
rier components in the principal solution are nearly tree of 
numerical dissipation. On the other hand, short-wavelength 
components may decay rapidly. 

(c) For a fixed e, Eq. (3. 1 8) implies that the principal solution 

is more diffusive for a smaller \v\. How to compensate this 

effect is a subject to be discussed in Section 7. 

(d) Equations (3.17) and (3.18) imply that, for all v with 

if < 1 and all 0 with -it < 8 =£ ft, we have 

()<E*,(e, v, tf)=£e, 0 < * (s, v, 0) < min{ 1 , 4e}, (3.19) 

which, according to Eq. (3.16), is equivalent to 

I - 8 < |G’ , r'| — E 1 — min{l,4c} ^ — 1 (3.20) 

As a result, by choosing e small enough, both |C:’| and |G’ UI | 
can be confined within an arbitrarily narrow range. As noted 
previously, the spurious part of a numerical solution generally 
is insignificantly small assuming a smooth initial condition. It 
does not contribute to accuracy and usually dissipates faster 
than the principal part. Thus, our primary concerns is how the 
principal part dissipates. From Eq.(3.20), one concludes that, 
for any e with 0 < e < 1, |C + 2| | will be bounded uniformly 
from below by a positive number 1 - s for all v with v < 1 
and all 0 with it < 8 ^ it. By choosing an e of proper 
magnitude, one can suppress artificial numerical oscillations 
without causing large diffusive errors for any combination of 
„ an d o. This fact contrasts sharply with what one expects from 
typical classical schemes which are usually very diffusive with 
respect to certain Mand 0, while not at all with respect to other 
r-and 0. As an example, we consider the Lax-Wendrott scheme 
(12. p. 101 1. Its amplification factor is of unit magnitude, for 
all 8 at v = 0, or v = 1. On the other hand, the amplification 
factor = 0 if ff = i and 9 = tt. 


In nonlinear flow solutions, e.g., shock-tube solutions to be 
discussed in Section 7, analogues of v are dependent on local 
velocity components. Thus, they may vary from one location 
to another. Also, at some neighborhood, the Fourier spectrum 
of the local solution may have peaks spread over a wide range 
of 8. Thus, for a numerical analogue of Eq. (2.22), a large 
variation in numerical diffusi vity with respect to (land v gener- 
ally means that numerical solutions obtained using its nonlinear 
extensions will suffer annihilations of sharply different degrees 
at different locations and different 8. Such selective annihila- 
tions may cause large distortions of numerical solutions [25]. 

This completes the discussion of stability and numerical 
dissipation. Other key subjects, i.e., consistency and the trunca- 
tion error, are discussed in Section 7 of [5]. 

In conclusion, the a-e scheme has been constructed to solve 
Eq. (2.22). It has the unique property that numerical dissipation 
can be controlled by a parameter e. Because neither characteris- 
tics -based techniques nor knowledge about the upwind direc- 
tion is used in the construction of the a-e scheme, as will be 
shown in the next section, it can be easily extended to model 
the Euler equations. 

4. THE EULER SOLVER 

We consider a dimensionless form of the 1 D unsteady Euler 
equations of a perfect gas. Let p, v, p. and y be the mass 
density, velocity, static pressure, and constant specific heat 
ratio, respectively. Let 

mi = p, u 2 = pv. w> = pHy “ 1 ) + (5)P y2 - (4 - 1 ) 


f = (y — 1 )«a + (!)(3 y)(ui)'lu\, (4.3) 

and 

f\ = yu 2 udu\ — (iKy — 1)(« 2 ) /(Mi) - . (4.4) 

Then the Euler equations can be expressed as 

+ = m = 1,2,3. (4.5) 

dt dx 

The integral form of Eq. (4.5) in space-time E 2 is 

i h„, ■ els = 0, m = 1 , 2, 3, (4.6) 

J ,S ( V > 

where K - (U ^ 2 < 3 < are lhe space-time mass, 

momentum, and energy current density vectors, respectively. 
As a preliminary, let 

= iff 8>h , m.k = 1,2.3, 


(4.7) 
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and let F be the matrix formed by f„ h . m, k = 1.2, 3. Let c 

be the sonic speed. Moreover, for any numbers ih a„, 

let diag(</,, a-., .... a„) denote the diagonal matrix with 

lh - a « being the diagonal elements on the first, second 

and nth rows, respectively. Then there exists a 3 X 3 matrix 
G such that 


Because h„, — (/,„ we also assume that 


h*(.r, t: j. n) = (/*(a. r, y. n), u*{. v. t\ j. n)), 
m = 1 , 2, 3. 


(4.15) 


G 'FG = diag(y, v - c,v + c), (4.8) 

where G 1 is the inverse of G. Note that v , c, F, G, and G" 1 
are functions of u„„ m = 1, 2, 3. These functions are given 
explicitly in [5|. 

Consider SEs of type I depicted in Fig. 2. For any (a, t) G 
SE( j, it). f, /). /„(.v. t ), and h,„(.r, t) are approximated by 

w *( v ’ r - .A «)./*(a, /; y. n), and h*(.v, /; y, «), respectively. They 
will be defined shortly. Let 


Note that, by their definitions: (i) ( fj " and (/„,)''. m = 1 , 2, 3, 
are functions of («,„)", tn = 1, 2, 3; (ii) ( /„,)", /» =1,2, 3, are 
functions of (h,,,)' 1 and («,„)", m = 1,2, 3; and (iii) (/„„); are 
functions of («,„); and («„„)", in = 1 , 2, 3. 

Moreover, we assume that, for any (a, t) G SE( y, n), u,„ = 
u*(x, t ; y, n) and/,, = /*(a, t\ j, n) satisfy Eq. (4.5); i.e., " 

<tu*(x,rj,n ) <■'!/’* (a, t‘,j,n) 

r b - = 0. 14 I At 


H,?(A, /; y, n) = („„,); + ( u „)J(x - a,) 

+ («„,)"(/ - f"). m= 1,2,3, 


According to Eqs. (4.9) and (4.14), Eq. (4.16) is equivalent to 


(4.9) 


( «*). j = -(/,,)"■ 


(4.17) 


where (i/„,)". («„„)". and (//„„)" are constants in SE(y, n). Obvi- 
ously, they can be considered as the numerical analogues of 
the values of hu m /hx. and du m /<H at (x : . /"), respectively. 

Let (/„)" and ( /„,;)" denote the values of/,, and / a , respec- 
tively. when ti m . m = 1. 2, 3, respectively, assume the values 
of (//,„)", in = I. 2, 3. Let 

3 

(./»,« );' = x (/m )?(«*, )", wt = 1 . 2. 3, (4. 1 0) 

A - 1 


and 


Because (/,„)" are functions of (u„y; and (u m )% Eq. (4. 17) im- 
plies that (m,„,)" are also functions of and («„„);. From this 
result and the facts stated following Eq. (4.15), one concludes 
that the only independent discrete variables needed to be solved 
in the current marching scheme are («,„)" and («,,„)/ 

From Eq. (4.16), one concludes that the generalization of 
the potential function /(a, r. y, „) introduced in Section 2 to 
the current solver are ifijx, r. j. n). m = 1, 2, 3, which satisfy 


bipjx, /; y, n) 

hi 


= /*( A, r; y, «) 


(4.18) 


3 and 

(,/m/ )/ ~ 2) (.///u )/'(wu )/*« fti — 1 , 2, 3. (4. 1 I ) 

I - ! 


Because 


y, aQ 

dr 


w*(.v, n). 


(4.19) 


d\* 


2 /. 


<>lh 

m. k 

dv 


Substituting Eqs. (4.9) and (4.14) into Eqs. (4.18) and (4.19), 

(4. i 2) and usin 8 (4. 1 7), one concludes that, up to an arbitrary con- 
stant, 


and 


i )t 




(4.13) 


(./mv); and ( /„„)" can be considered as the numerical analogues 
of the values of <)f m U)x and <{f„J <H at (a,, t"), respectively. As a 
result, we assume that 


iA„,(a, r; y, n) = (/,)”(/ - /") - (x - a,) 

+ (Mf„,y;u - n 2 - (bu< m j;(x - x , )- (4.20) 

+ (/«,)"( a - a ,)( r - /"). 

By using an argument similar to that leading to Eq. (2.36), one 
concludes that 


/*(a, t: j, in = ( f„y; + - Xj ) 

+ (/„,)"(/ - t"), m = 1,2,3. 


(4.14) 


j, h * ' d* ~ <h(x\ t’ J, n) - ifijx, r; / n). (4.21) 

Here I is a simple curve joining (a, t) and (a', /'), and lying 
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entirely within SE(./, n). We also assume that (is points to the 
right of I' if one moves forward from (x, t ) to t'). 

As in the u-e scheme, we assume that the flux of h,t is 
conserved over CE(/\ n ), i.e., 


marching scheme presented in [3] is formed by Eqs. (4.24) and 
(4.28) with s = i 

To construct a larger class of generalizations to Eq. (3.1 Oh 
for all (./. n) E il, let 


SdMijjt)) 


(Is = 0. 


Combining Eqs. (4.21) and (4.22), one has 

<//„,(.v, - A.v/2, t": j. n) - tfrjx , + Ax/2, j. n) 

+ >pjx, + Aa/2 , /" l,: ; j - i « - 5) 

- </av, l/: + Af/2;y-i«-i) 

+ i hiXjon.-t" " i + Ati2:j + ln-\) 

- <lrJ.x ,. U 2 - At/2, r l,: ; j + ill- h) = o. 


(4.22) 




m= 1,2,3. (4.29) 


Let (ej’\ /n = 1, 2, 3, be parameters that can be functions of 
m — 1 , 2. 3. There can be many choices ot these functions. 
Let ( £v)“ be the value of the On, A) -element of the matrix 6 
when m - 1 . 2 , 3 , respectively, assume the values of 

(uj-;, = 1, 2, 3. Similarly, let ($J); be the value of the 

(4 ? 3 ) (in- lfc)-elemcnt of the matrix G 1 when m = L 2, 3. respec- 
tively, assume the values of 00 / • m = 1 • -• - 4 - ^ et 


= '"•* =I - 2 - 3 - (430) 


Substitution of Eq. (4.20) into Eq. (4.23) yields 


00" = iioo; Hi + (« Jr |V 2- + (*-)?• i/i - in: I- < 4 - 24) 


where, for all ( 7 , m) E 11, 


+ «= 1.2.3. 


(4.25) 


Equation (4.24) forms the Hrst half of the current marching 
scheme. The second half which solves (m„„)” will come from 
a generalization of Eq. (3.10). 

For all (j. n) G il, let 


Then Eq. (3.10) can be generalized as 


(M mi ); = [<«:„)"- 1/2 - 04)./ i :]/A v 

3 

+ ^ (2(0 )" — <5„ l t|((/uu) ) . 


(4.31) 


where in = 1 , 2 , 3 , and <?,„* is the kronecker-delta symbol. 

Consider the special case in which, tor all ( ./. ») G Q- 
(£,)" = (e : ); = (e,)?. Let (e,„)>- («)'/■ '» = >• 2 - 3 - Then 
(e„,*)" = (e)''5 mi , and thus Eq. (4.31 ) is reduced to 

(Mm,)" = I 04 )/ , 1/2 - (M,„)" i ;|/Av ^ ^2) 

+ 12 (e)" - 1 \Uhl,„Jn m = l - 2 - 3 - 


Ulll„n)'/ — 2|(M„,.v)'/m I / 2" + (M„„ )/-l/2 1 ^ 

- I (Mm)", Hi - (M„,)rml/A.r 

and 

(«i)J, i /2 = (uJ'liHi + (A//2 )(«,„/);* (4-27) 

for m = 1, 2, 3. Because Eqs. (4.26) and (4.27) are the general- 
izations of Eqs. (3.2) and (3.11). respectively, a natural general- 
ization of Eq. (3.10) is 


Note that Eq. (4.32) reduces to Eq. (4.28) if (£)7 e f° r a11 
( /, w) E 11. 

Recall that both y and <■ are functions of u m m = 1, 2, 3. 
For all SE(y, n), let 0" and c”, respectively, denote the values 
of it and c when u,„, m =1,2. 3, respectively, assume the 
values of (m„,)". m = 1,2, 3. It is shown in |5] that the marching 
scheme formed by Eqs. (4.24) and (4.31) is stable if, for all 
( j. n) G 12, 


where 


(m„j; = [(«:);, i/2 - (m,'„)/-i/-|/a.v (4 28) 

+ (2s - 1 )(r/«,„,)", in = 1,2,3, 

where e is a parameter independent of numerical variables. 
Note that the last term in Eq. (4.28) vanishes if e = i The 


We conclude this section by introducing some possible modi- 
fications to the above solver. Note that 04)". by its definition. 
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n + 1 
n + 1/2 
n 

n - 1/2 
n - 1 


(J.n) 


(J “ 1/2, n - 1/2) (j + 1/2, n - 1/2) 


but not between (y, n ) and (y — 4, /?), one would expect that 
> |(«i,u-)y|. Moreover, because (y, n) and (j - i n) 
are on the same side of the discontinuity while (y, n) and 
(j + i «) are on the opposite sides, (u mx )” should be a weighted 
average of (u n ^)? and biased toward the one with the 

smaller magnitude. 

As a result of the above considerations, (w; m )' r can be re- 
placed by 


UC; )" = W lt ((u mx _ )'', (u mx , );; a), m = | , 2, 3, (4.38) 

Here a is an adjustable constant and the function W„ is defined 
by (i) U/,(0, 0, a) = 0 and (ii) 


(4.39) 

(kl + \x. I > 0), 


FIG, 5. The mesh and CEs of Ihe u-f scheme, (a) The relative positions 
of CEs and mesh points, (b) CE{ j. n). 


represents a finite-difference approximation of u m at (j ± n). 

As a result, 


(w«i); 1/2UA.V, m — 1,2,3, (4.35) 

respectively, are the central-difference approximations for 
hu m f h.\, m — 1, 2, 3, at (y, n). Note that («(,„)•' is the first term 
on the right side of each of Eqs. (4.28), (4.3 1 ), and (4.32). The 
above central-difference approximation is valid as long as no 
discontinuity ot u m (or its derivatives) occurs between (j — 5, 
n) and (y -f n) (see Fig. 5). In the following discussion, 
we develop alternates which are valid even in the presence 
of discontinuity. 

Let 


)/ - 1 /2 (Um)j 

' “ IIJl ' 


1,2,3, (4.36) 


where (u m y; can be obtained from Eq. (4.24). Because 
( l *m)/ )/, and are the numerical analogues of 

14 m at (J ~ 2, «), (y, n) and (y + 4, n) % respectively, 
(u mx )■ and (« mit )■' are tw-o numerical analogues of the value of 
dujte at (y, n) with one being evaluated from the left and 
another from the right. Note that 

( u nx )j ~ 2 f( Umx )j T ( Unix +)/']■ (4.37) 


where x- and x- are any two real variables. Note that W t) (x , 
-v,; a) = (x + xJ/2: i.e., («-); = («<j;, if a = 0 or |/| = 
Also the expression on the right side of Eq. (4.39) repre- 
sents a weighted average of x and x. with the weight factors 
|.v.|‘V(|x + |" + \x | a ) and |x-| tr /(|x t | w + |.v |"). For a > 0, this 
average is biased toward the one among a\ and x with the 
smaller magnitude . For the same value of [\\\ and \x |, the 
bias increases as a increases. Thus, we should always choose 
a > 0. 

Note that the special weighted averages Wjx , x 4 ; 1) and 
VT^U-, x i ; 2) are used in the slope-limiters proposed by van 
Leer (28) and van Albada 1 29] , respectively. 

The above modification, i.e., («;,j; replaced by («;;;;);, is first 
given in [3 J. It is shown in [3 J and also Section 7 of the current 
paper that it is an efficient tool to suppress overshoots and/or 
numerical oscillations near a discontinuity. Moreover, because 
7g u + )j are constructed using only the data associated with the 
mesh points (j - J, « - 4) and ( j + i n - 4), the effect of 
this modification is highly local; i.e., it generally will not cause 
the smearing of shock discontinuities. 

However, there may be a price to pay for the above modifica- 
tion. Because a fractional power is costly to evaluate, so is 
lV f> (x , a + ; ex) if at is not an integer. Moreover, because the bias 
of this weighted average increases with a, a situation may arise 
such that the use of an a with |a| < 1 may be desirable. To 
obtain a computationally efficient weighted average of arbitrary 
small bias, let 


W(.x , .v. ; a. 13) = ( 1 - (i)Wjx , A . ; 0) 
+ PW'Xx.. x . ; a). 


(4.40) 


In case a discontinuity occurs between (y, n) and (y + £, n) 


where ft ^ 0 is an adjustable weight factor, and a generally 
is an integer. Because W„( x., x*; 0) is the simple average of 
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jt and jc+, Eq. (4.40) defines a linear weighted average of this 
simple average and the nonlinear weighted average defined in 
Eq. (4.39). Obviously, W(x-, x^; a, /3) = (i)(x_ + x+) if x . = 
x+. Furthermore, because 


tom ffjn _ 
dt dx 2 


m 


1,2.3, 


(5.4) 


The integral form of Eq. (5.4) in space-time E 2 is Eq. (4.6) with 


W f ,(x_,x + ; ~ot) = 


| x-|"x+ + l A \ X- 
|x,|“ + \x~\ a 

(W + k i > 0). 


(4.41) 


h,„ = (/„ - afjax, u m ), m = 1 , 2, 3, 

As a preliminary, let 


(5.5) 


alternatively, W(x , x f ; a, fi) can also be expressed as 


f m .k = m, k = 1,2, 3, (5.6) 


Wu, x+; a, p) = | ^>( A -» A ’ a ) 

' (4.42) 

+ (L=i?j w -< x 

The application of the more general modification, i.e., OCO/ is 
replaced by 


(«y; = W((u mx Yf (iW,)?; a, P), 2, 3, (4.43) 

will be demonstrated in Section 7. 

Finally, note that W( x_, x+ ; a, /3) can be further generalized 
by a linear weighted average of several Wf x , x> \ a) with 
different values of a. 

5. THE NAVIER-STOKES SOLVER 

We consider a dimensionless form of the ID unsteady Na- 
vier-Stokes equations of a perfect gas [12, pp. 191 — 1931. (Note: 
the expressions on the right sides of the last three equations in 
Eq. (5-47) of [12] have incorrect signs in the earlier versions. 
The conduction heat-flux vector should be proportional to the 
negative of the gradient of temperature.) These equations are 
extensions of the Euler equations defined in Section 4. Thus, 
unless specified otherwise, the symbols, definitions, and equa- 
tions given there will be used in this section. 

Let Re L and Pr denote the Reynolds number and Prandtl 
number, respectively. They are assumed to be nonnegative 
constants. Let 


/'=-<>. 

~ def 4 a 2 

3 Re, Mj ’ 


and 


2? UCI 

" ~ 3 Re, 



Re, Pr 


U[ 


( H ;) 2 

2 ( w ,) 2 


(5.1) 

(5.2) 


(5.3) 


and 


dcf 4 y 

~ 3 Re,,’ Tl ~ Re, Pr 


def 

Ty = J 2 — T, . 


(5.7) 


Let F denote the 3 X 3 matrix formed by f , u , w, k = 1, 2, 3. 
Then Eqs. (5. 1 )— (5.3) imply that 


F = 


I 


0 

TlM. 

( M ,) 3 


1 (Mi) 2 Uy 

T2 (M,) 2 


0 ()\ 

— 0 

U\ 

Tilt 2 Ti 

(M,) 2 1* if 


(5.8) 


Again we consider SEs of type I depicted in Fig. 2. For any 
(x, t) E SE(/ nf M m (x, tffjx, t)Jjx , /), and h,„(x, /), respec- 
tively, are approximated by w*(x, t\ 7 , /?), /*(x, /; 7 , w), ?*(x, 
r; 7 , n), and h*(x, f; 7 , a?); m*(x, f; 7 , /i) and /£(x, /; 7 , /?), 
respectively, are defined in Eqs. (4.9) and (4.14);/*(x, t\ 7 , n) 
and h*(x, t\ 7 , n) will be defined immediately. 

Both /,, and f mk are functions of w„„ m = 1, 2, 3. Let (/„)" 
and (/„.*)", respectively, denote the values of/„ and f mk when 
m„„ m = 1 , 2 , 3 , respectively, assume the values of (u m )'f m = 
1, 2, 3. Let 

(.L)/ = 2 (,/L )"(«u)". iw = 1 , 2, 3, (5.9) 

fc- I 

and 


(/»>; = 2 (/-.*)?(«*,)?. = 1 ■ 2. 3. (5.1 0) 

Jt^l 

Using an argument similar to that leading to Eq. (4.14), we 
assume that 


/*(.r, t: j, n) = (/„,)," + - x } ) 

+■ (t ~ t"), m= 1,2.3. 


(5.11) 


Then, the Navier-Stokes euations can be expressed as 


As a result of Eq. (5.5), we also assume that 
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h*(.V, t\ /', n) 

( rV*(.v, y, n) \ 

= — ,i<2U,/;y,/i)J, (5.12) 

m = 1.2,3, 

Also, we assume that, for any (a, /) E SE( /, n ), it m — 
//*(v. K y, n), = /*(a\ /; y, /?), and X, = /*(a. /; y, u) 

satisfy Hq. (5.4), i.e.. 


du*(.v, /: y. /?) 


(U 


d 

4- — 

rfv 


(5.13) 


/*(a\ f; y, n) - 


/; y, /i) 


d.V 


= 0. 


The above condition again leads to Eq. (4.17). Thus, the diffu- 
sion term in Eq. (5.4) has no impact on how m*(a\ /; y. /?) varies 
with time within SE(y, n). This same fact was observed in 
Section 2. The reason behind it and its significance were also 
discussed there. As a result of Eq.(4. 17), and other definitions 
given earlier in this section, one can conclude that the only 
independent discrete variables needed to be solved in the current 
solver, as in the Euler solver described in Section 4, are also 
(u,„y; and 

A comparison between Eqs. (4. 16) and (5. 1 3) reveals that, for 
the current solver, Eqs. (4. 1 8) and (4.19) should be replaced by 


(f m )" in Eq. (4.20) is replaced by (/,,,)• in Eq. (5.16). Obviously, 
Eq. (4.21) is still valid for the current solver. Because i 
f; y, n) is independent of (/„,)■ and (f M )% Eq. (4.2 1 ) implies that 
the last two parameters are irrelevant in flux evaluation. More- 
over, because the current solver will be constructed using only 
flux-balance conditions, these parameters are also irrelevant in 
the following construction. 

For all (y, n) E II, we assume that 

i h*-</s = (). (5.18) 


With the aid of Eqs. (5. 16) and (4.21 ), Eq. (5.18) implies that, 
for all (y, n) E II, 

(«»,)" “ ± ^ !(»„„ Yj'-'i: + (w„j;| 

± W,„)'F Ms - (U);i (5.19) 

( A t U 

± 4A7 |( /”" ) '' |l,j! + = °- 

Adding the two equations given in Eq. (5.19) results in 

(«„,)" = *l(«ori« + («»)?. Ms + (•<■»,)'/ Ms - (Sjj-Hsl (5.20) 

where, for all ( j, n) G li. 


,h fa,(x.r;j,n) . i)f*(x, t\ j. n) 

■ =J£(x. t: j, n) — (5.14) 


lit 


dx 


and 


, det 4 A 4 1 

(.s,„); = —(«„„); + —(./,„); 

+ m= 1,2,3. 


(5.21) 


respectively. Note that Eqs. (5.15) and (4.19) are identical. 
According to Eq. (5.1 1), the second term on the right side of 
Eq. (5.14) is simply the constant — (/, u )". Thus, for the current 
solver. Eq. (4.20) should be replaced by 


Equations (5.20) and (5.2 1 ) are the current counterparts of Eqs. 
(4.24) and (4.25), respectively. By using Eq. (5.20), (u m ) n } can 
be solved explicitly in terms of discrete variables at the next 
lower time level. 

By substruction of the two equations given in Eq. (5.19) and 
using Eq. (5.17), one has 


t/Uv. K y, n) = ( - t fl ) - (mJJ’U - a- ) 

+ (*)(./;, ,,)"(/ - l"Y - (3)(H„„)"(.v - Xj) : (5.16) 

+ (/,„)!'(.'■ - xj)(t - 1"). 


At 




(.ly; 


(A/) : 

4Aa 

+ f~ x <.L)" = </>,,, )". m = 1,2.3, 


(5.22) 


where 


where, for all (j. n) G (i, and m = 1, 2, 3, 


(fjj = (/»)" - (L)"- (5.17) 

The only difference between Eqs. (4.20) and (5.16) is that 


(bj'j = (/„)" + ?| («,„)", Hi 

~ (UmYj Mi ~ (Oyt Ms — (•V»)J Mi\- 


(5.23) 
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Note that (/„)•, m = 1, 2, 3, are functions of (mJ>\ m = 1,2, 
3, and the latter can be evaluated by using Eq. (5.21). Thus, 
{b, n )% m = 1 , 2, 3, can also be evaluated in terms of the variables 
at the (n - i)th time level. 

To proceed, note that Eqs. (4. 10), (4.1 1 ) and (4. 17) imply that 

.3 3 

(/*)? = -22 (5.24) 

it - ! /“ I 


where q*(/i, 0) is a 2 X 1 column matrix. Substituting Eq. 
(6.1) into Eq. (2.18), one obtains 

q*(n + 1, 0) = \Q( v, £ fl)! : q*(«, 0), (6.2) 

where 

Q( v , £ 0) = e~' m Q + + e"" 2 0 . (6.3) 


Moreover, for all (y, /t) E 11, let 


(«„'„ )" = ^ («,,,);. »» = i . 2 , 3, 

4 


(/,;,„ 

yt ^ ( f yi 

k), ~ K yhiUJi' 

m, k = 1 

, 2, 3. 

At 




del 4Af 

, /», k = 1 

1,2, 3, 


k) > _ (A.V )- Umkh 


and 


(cu)] = A* + (£,*)? 


2 '”>* = >.2.3. 


(5.25) 

(5.26) 


According to Eq. (6.2), the amplification matrix is the square 
of the matrix £7(>, £, 7). Substituting Eqs. (2.16) and (2.17) 
into Eq. (6.3), one has 


GO',*,#) 


A/ii 7i: 
V ^/21 722 


where 

(5.27) 

cy i , = cos(tf/2) - />sin(0/2) 
cj ] 2 = -i(l - “ f) sin(0/2) 

/(I - v 2 ) sin(0/2) 

(5 28) a,, = - — f[cos(0/2) + iv sin(0/2)). (6.4) 

1 — IT + t; 


Let 


With the aid of Eqs. (5.9) and (5.24)— (5.27). Eq. (5.22) can be 

reexpressed as V(v, £ 0) = £cos(0/2) - />( 1 - v 2 ) sin(0/2). (6.5) 


2 («».*)/(«*,)" = (£»)". w =1,2. 3. (5.29) 

k - l 

Because (/„, ,)" and (/,;,()", m, k = 1,2, 3, are all functions of 
(«,„)", tu - 1,2, 3, so are (a mJ )" m, ( = I, 2, 3. Thus, («„*)" 
can also be evaluated in terms of the variables at the ( n - s)th 
time level. It follows that, for each (j. n) £ li, Eq. (5.29) 
represents a system of three linear equations for three unknowns 
(«,;,);, m — 1, 2, 3. These unknowns (and thus («„„)", m = 1, 
2, 3, through Eq. (5.25)) can be solved easily by a matrix 
inversion. Equations (5.20) and (5.29) form the current 
marching scheme. 


Then the eigenvalues of Q( v , £ 0) are 
v-(v, £ 0) 

« y(v, £ 0) ± \ !>)(<'•£ 0W + (1 ~ ^) 2 ~ ? ( 6 ' 6) 

1 - p 2 + £ 

Thus the amplification factors GV’and G_’ of the «-/x scheme 
are given by 

CL" = [o-.(^, £0)p. (6.7) 

Note that 


6. STABILITY ANALYSIS 


G 1 , 1 ’— » I, 


I - v 2 - g 
!-(/- + £ 


as0-*O (6.8) 


The stability of the a-yr and a-e schemes will be studied 
using the von Neumann analysis. For all (/, /i) £ fi, let 

q( j, n) = q*(n, 0)e' J " (/ = V— 1, —tt< 0 ^ n), (6. 1 ) 


if 1 — jr > 0. Because the amplification factor of a plane- 
wave solution to Eq. (2.1) approaches 1 as 0 — » 0, G'J 1 and 
G"’ are referred to as the principal and the spurious amplifica- 
tion factors, respectively. Moreover, Eqs. (6.5)— (6.7) imply that 
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Gl" = |y-|—|[|cos(0/2) - /Vsin(0/2) 


± V||cos(0/2) - f>sin(0/2)] 2 + 1 - I 2 ] 


(6.9) 


M(e , v , 6) 
/cos(0/2) 
~\ i(l 


ivs\n(0/2) 
e)sin(0/2) (2e 


-/(l - i^ 2 ) sin(0/2) 
l)cos(0/2) - fVsin(0/2)> 
(6.15) 


if 1 — ir ^ 0, and f = £/(l — v 1 ). Similarity between Eqs. 
(6.9) and (2.21) was noted in Section 2. 

In 1 1 ], the stability of the a-jji scheme is studied using a 
rigorous discrete Fourier analysis. The von Neumann stability 
analysis can be considered as a limiting case of the discreate 
Fourier analysis. By using Eqs. (4.29) and (4.30) in [1], one 
can infer that the a-fi scheme is stable if and only if, for all 6 
with 77 < 0 < 7T, 

max{|GV h |, |G 4,1 |} ^ 1 if Q(v, £ 0) is nondefective (6.10) 
and 


The eigenvalues A~(e, v , 6) of M(e , v , 0) were given in Eq. 
(3.13). The principal amplification factor G r p and the spurious 
amplification factor G_ ) of the as scheme were given in Eq. 

(3.12). Note that 

G <2> — ► I, G 2 _^2e— 1 as0— ►() (6.16) 

if Eq. (3.14) is assumed. Moreover, from Eqs. (6. 1 0) and (6. 1 1 ), 
one infers that the as scheme is stable if and only if, for all 
0 with —7i < 0 < n, 

max{|G (2) |, |G' 2i |} < 1 if M(e. v , 0) is nondefective (6. 17) 


|Gl n | < 1 if Q( v , £ 6) is defective. (6. 1 1 ) 

Note that G { V = G lh if Q{v , £ 0) is defective |30, p. 353]. 
Assuming £ > 0 and 1 v' + £ ^ 0 (the latter is a basic 
assumption of Eq. (2.14)), it is proved in [1] that the current 
scheme is stable if and only if v~ < 1. 

Let (1 - irj 1 7^ £ 2 such that both Eqs. (2.14) and (2.24) are 
valid. Combining Eqs. (6.5)— (6.7), one has 


GVG' 



(6.12) 


and 


|G f r1 < 1 ifA/(c, i\ 6) is defective. (6.18) 

Equation (3.13) implies that 

|A + (e, v, 0)| |A_(e, »/,0)| = \2e - l|. (6.19) 

By using Eqs. (3.12) and (6. 17) — (6. 19), one concludes that 
stability requires that \2e - l| < 1, i.e., 0 < e < L Thus the 
first part of Eq. (3. 14) is necessary for stability. Equation (3. 1 3) 
also implies that 


Because the amplification factors of the backward-marching 
scheme are (G ( + n )' 1 and (G m ) ' stability of both Eqs. (2.14) 
and (2.24) requires that | G | = |G. h | = 1. According to Eq. 

(6.12), the last condition cannot be met if /jl > 0 and \P ^ 1. 
This result was used in a discussion given in Section 2. 

Next we study the stability of the as scheme. By substituting 
Eq. (6.1) into Eq. (3.7), one has 

q *{n + 1, 0) = \M{f„ v, 0)| 2 q(n 0), (6.13) 

where 

M(f„ v, 0) = e ~ m M , 4- e‘ m M . (6. 14) 


A.(e, v, tt) = —iv± V(l - e)(l - v 1 ). (6.20) 

Thus, 

max{|A,(e, i\ 7f)|, |A_(e, v , 7r)|} >1 if iP > 1; e < 1. (6.21) 

The first part of Eq. (3.14) coupled with Eqs. (6.17), (6.18), 
and (6.21) implies that iP < 1 is necessary for stability. Because 
the case v 1 = 1 is ruled out by the basic assumption l — 
iP # 0 of Eq. (3.6), the second part of Eq. (3.14) is also 
necessary for stability. The proof that Eq. (3.14) is also suffi- 
cient for stability will be given later in this section. 

To prove Eqs. (3.17) and (3.18), note that Eq. (3.16) im- 
plies that 


According to Eq. (6.13), the amplification matrix of the as AT( e * v -> 0) “ e (x’ ■+■ (6.22) 

scheme is the square of the matrix M(e, v , 0). Substituting Eqs. 

(3.8) and (3.9) into Eq. (6.14), one has 


where 
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X' = (1 - f 2 ) sin 2 (0/2) + 2(1- e) cos : (0/2) (6.23) 

and 

X" = 2 cost 0/2) 

X V(1 — e)|(l — e) cos 2 (0/2) + (1 — f 2 ) sin _ (0/2)[. 

(6.24) 

With the aid of Eq. (3.14) and -jt < 0 <= jr. , Eqs. (6.23) and 


(6.24) imply that 

X' = x" = 0 if e = 1 ; 0 = 0, (6.25) 

f =0, if 8 = 1; 0 = 0; 

y ' \ (6.26) 
{>(), if f, 1 ; or 0 + 0, 

X" >2(1 - e) cos 2 ( 0/2) > 0, (6.27) 

x' — x" — ( 1 — f 2 ) sin 2 (0/2). (6.28) 

(V ~x"M + x") = (x')- -(/)-’ 

= (1 - f 2 ) 2 sin J (0/2). (6.29) 


For the case e = 1 and — 0, Eqs. (3.17) and (3. IK) follow 
immediately from Eqs. (6.22) and (6.25). 1 hus, in the following 
proof of Eqs. (3.17) and (3.18), we assume that 

e 1 or 0 # 0. (6.30) 

Combining Eqs. (6.26), (6.27), and (6.30), one concludes that 

X + X" > (6.31 ) 

Equations (6.29) and (6.31) imply that 

X 9 - X* ^ o. (6.32) 

Equation (3.18) now follows from Eqs. (3.14), (6.22), (6.28), 
and (6.32). The validity of the first inequality sign in Eq. (3.17) 
follows from Eq. (3.18) and the tact that 8(1 “ s) — 0 it 0 — 
£ < 1. The validity of the second inequality sign follows from 
the fact that 

* (e, f. 0) - *-(e, f. 0) = lex" (6 33) 

> 4e( 1 - e)cos 2 (0/2). 

Equation (6.33) is a simple result of Eqs. (6.22) and (6.27). To 
establish the validity of the last inequality sign in Eq. (3.17), 
note that 


X (s, f, 0) = e(*' + x") = e \ 2 X' ~ _ X")\ - 2r X' 

= 2e|(1 — F 2 )sin : (0/2) (6 34) 

+ 2(1 - e)cos 2 ( 0/2)| 

< max{2e( 1 - f 2 ), 4 e( 1 - e)} ^ 4 b, 

where Eqs. (6.22), (6.32), (6.23), and (3.14) have been used 
Moreover, because |G‘ 2, | s 0, Eq. (3.16) implies that 

X (e, f, 0) — 1 • (6.35) 

The validity of the last inequality sign in Eq. (3. 17) now follows 
from Eqs. (6.34) and (6.35). Q.E.D 

Next we shall prove that Eq. (3.14) is also sufficient for 
stability. Note that, as a result of Eqs. (3.17) and (3.18), 0 - 
^.(B, f, 0), and thus |G'3’| ^ 1, for all e, f, and 0 satisfying 
Eq. (3.14) and — 7 r < 0 ^ jt. As a result, Eq. (6.17) is always 
satisfied. To complete the proof, we need only show that 
Eq.(6. 1 8) is also satsfied. To proceed, note that G 1 , 2 ' = G' ' it 
A7 (e, f, 0) is defective. From Eqs. (3. 1 2)— (3. 14), one also 
concludes that s = 1 is necessary if G 1 . 2 ' = G' 2 '. Moreover. Eq. 
(6.15) implies that M{ 1. f, 0) is the identity matrix. Thus, one 
concludes that e — 1 and ^ / 0 are necessary if A/(e, f, 0) is 
defective. Because (i) 

G< 21 = [cos(0/2) - if sin(0/2)] 2 if e = 1 (6.36) 

and (ii) 

|[cos(0/2) - /Fsin(0/2)] 2 | <1 if f 2 < 1 ; 0 ^ 0, (6.37) 

one arrives at the conclusion that Eq. (6.18) is also satisfied. 

Q.E.D 

7. NUMERICAL RESULTS 

In [ 1 ], numerical solutions of Eq. (2.1) generated by the 
MacCormack [12, p.102], the Leapfrog/DuFort-Frankel, and 
the a-fjL schemes are compared with the corresponding analyti- 
cal solutions for different values of physical coefficients, mesh 
parameters and total marching times. These comparisons show 
that the a-fi scheme is far superior to the Leapfrog/DuFort- 
Frankel scheme in accuracy and has a substantial advantage 
over the MacCormack scheme in both accuracy and stability. 

In this section, accuracy of both the Euler and the Navier- 
Stokes solvers will be evaluated numerically using a shock tube 
problem suggested by Sod [31). Because the u-e scheme may 
be considered as a special case of the Euler solver, no separate 
numerical evaluation for the a-e scheme will be given. 

Let the specific heat ratio y = 1.4. At / = 0, let (i) 
(p, v, p) = (I, 0, 1), i.e., (Mi, ii :■ U\) = (1.0, 2.5) it x < 0, 
and (ii) (p, u, p ) = (0.125, 0, 0.1), i.e., (iq, M;, iq) — (0.125, 
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-0.4 -0.2 0.0 0.2 0.4 0.6 



0.4 - 0.2 0.0 0.2 0.4 0.6 


x 

FIG. 6. The Kuler solution (e = 0.5, a - 0. Af = 0.(X)4, CFL = 0.88). 


assumptions imply that the computation domain can be limited 
to hi £ 50. 

In the initial evaluation, we consider the Euler marching 
scheme defined by Eqs. (4.24) and (4.28). Numerical results 
(triangles) obtained assuming At = 0.004 and e = j are com- 
pared with the exact solutions (solid lines) in Fig. 6. Because 
each marching step advances the solution from Mo t 4- At/2, 
these results at t = 0.2 are obtained after 100 steps. Also it 
can be estimated that CFL = 0.88, where CFL is defined to 
be the maximum value of (|t/| + |o|)Af/Ax Thus the numerical 
calculation is carried out within the stability limits given by 
Eq. (4.33). Note that the agreements between the numerical 
results and the exact solutions are excellent. Particularly, the 
shock discontinuity is resolved almost within one mesh interval, 
and the contact discontinuity is resolved in four mesh intervals. 
Also, there are only slight numerical overshoots and/or oscilla- 
tions near these discontinuties. 

According to the discussions given in Sections 3, 4, and 6. 
the Euler solver behaves like the Leapfrog scheme, if e = (), 
and like the Lax scheme, i I & ~ 1. The former is free from 
numerical dissipation while the latter is highly diffusive. The 
current scheme with e = 5 can he considered as a scheme 
midway between the above two celebrated schemes. 

Moreover, the last term on the right side of Eq. (4.28) van- 
ishes if e = i The remaining term is simply a central-difference 
approximation for («„)". 


0. 0.25) if x > (). For all ( /. n) £ (1, let x t = jAx, and t" = 
nAt. Then (i) 

((«i>". («:)". («,))’) 

f( 1.0. 2.5) if 7 =-i -I,...; (7.1) 

1(0.125,0,0.25), if/ = i| 

and (ii) («„,)}’ = 0. j = ±i ±f, .... for m = 1,2, 3. Hereafter, 
we assume that n > 0. 

The above initial conditions coupled with several equations 
given in Sections 4 and 5 imply that, for both the Euler and 
the Navier-Stokes solvers, («„,)" is a constant and (u n ,)" = 0 
in two separate regions that are defined by j < -(« + i) and 
j 2: (n + 5), respectively. Thus, one needs to evaluate the above 
variables only if |y| < (n + 5). 

Without exception. Ax = 0.01 is assumed in this section. 
Also, all numerical results will be compared with the exact 
weak solution at t = 0.2. Because, at t = 0.2, the effect of the 
initial discontinuity at t = 0 is far from reaching the spatial 
regions defined by x > 0.5 and x < -0.5, respectively, numeri- 
cal computations, unless specified otherwise, will be simplified 
by assuming that, for all it with t" < 0.2, (i) 


((«,)", («.);', (Mv)") = 


(1, 0,2.5) 


if x < -0.5; 


(0.125,0,0.25) if x > 0.5, 



and (ii) (»„,,)" - 0 if |.v,| > 0.5. Because A.v = 0.01, the above 


FIG. 7. The Euler solution (e = 0.5. a = 1,4/ = 0.004. CFL = 0.SX). 
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X 


FIG. 9. The Euler solution (e = 0.7, « = 0, A r - 0.004, CFL - 0.88). 

FIG. 8. The Euler solution (e = 0.1, ct — 0, A/ = 0.004, CFL — 0.88). 


Let Eq. (4.28) be modified with (u^)* being replaced by 
(u^)j (see Eqs. (4.35) and (4.38)). Again assuming that A / = 
0.004 and e = J, the numerical results obtained with a = 1 are 
given in Fig. 7. The results obtained with ol 2, and ot 3 
are almost identical to those shown in Fig. 7 [51- The effective- 
ness of the above modification as a tool to surpress numerical 
wiggles near discontinuities is apparent. It was explained in 
Section 4 why this modification does not cause the smearing 
of shock discontinuities. Furthermore, the modification has no 
discernable effect on the smooth part of the solution. Because 
= (u l mx )] if a = 0, in the following discussion, it should 
be understood that the above modification is turned off if ot — 0. 

Note that the results shown in Figs. 6 and 7 can be generated 
using the sample program listed at the end of the present paper. 
It is coded assuming e = 0.5 and without imposing the condi- 
tions given in Eq. (7.2). The parameter a is represented by ia 
in the code. 

Let a = 0 and Af = 0.004. The numerical results obtained 
with e = 0.1, and e = 0.7, respectively, are given in Figs. 8 
and 9. Note that the case with e = 0.5 are given in Fig. 6. 
p or e = o. 1, because the scheme has very small numerical 
dissipation, pronounced wiggles appear in large regions near 
discontinuities. However, because of the same reason, the 
smooth part of the solution is highly accurate. The results 
shown in Figs. 6, 8, and 9, and other results obtained with 
e = 0.3 and e = 0.9 [51 are consistent with the theoretical 



-0.4 -0.2 0.0 0.2 0.4 0.6 



-0.4 -0.2 0.0 0.2 0.4 0.6 



0.4 -0.2 Q.C 0.2 


x 

L 10. The Euler solution (e = 0.5, a - 0, Af = 0.0004, CFL = 0.088). 
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X 


FIG. 11. The Euler solution (e = 0.5. « = 1, A/ - 0.0004, CFL = 0.088). 


prediction that the Euler solver becomes progressively diffusive 
as the value of e increases from 0 to 1. 

The above numerical results are all generated assuming 
At = 0.004. The numerical results shown in Figs. 10 and 1 1 
are generated with It = 0.0004 (i.e., CFL = 0.088). Note that 
now it takes 1000 marching steps to advance the solution to 
t = 0.2. Other defining conditions for these figures are identical 
to those tor Figs. 6 and 7, respectively. A glance over Figs. 6, 
7, 10, and 1 1 reveals that the current Euler solver is more 
diffusive at a smaller CFL. Note that, by considering the trunca- 
tion error, it was shown in |5| that, for constant e and Ax, the 
a-E scheme becomes more diffusive as A t decreases. A similar 
conclusion can also be reached by studying the amplification 
factors given in Eqs. (3.12) and (3.13). Because the Euler solver 
is a straightforward extension of the a-E scheme, one would 
expect that the former also behaves similarly. 

Also, as the value of CFL decreases, the diffusive effect of 
replacing a = 0 with a = 1 generally becomes more dis- 
cernahle. In other words, numerical dissipation introduced by 
replacing a = 0 with a > 0, is greater when CFL is small. 

To modify the above Euler solver such that it can compensate 
tor the observed effect of increasing numerical dissipation as 
Ar decreases, in the following discussions, we shall consider 
the more general marching scheme defined by Eqs. (4.24) and 
(4.32). The parameter (£)'' in Eq. (4.32) will be dependent on 
the mesh position (y, n) and the ratio At/ lx. Moreover, the 


term << v ); (see Eq. (4.35)) in Eq. (4.32) will be replaced by 
which is defined in Eq. (4.43). The weight factor £ will 
also be dependent on (y, n) and At/ Ax. 

To proceed, let 

v . del 

£(x) = x exp( 1 - .v), 0 < x < 1 . (7.3) 

Because £ is an increasing function within its domain, we have 
£U) <£(!)= I, 0<.r<l. (7.4) 

For all ( j, n) E if, let 

(£))" = bZUV^J!) (7.5) 

and 

(<,); = wt(u„ n y;,(u m + );: a, V(P, „„)"), (7.6) 

where ( tv,,)" is defined in Eq. (4.34), and b and a are constants 
that do not vary from one mesh point to another. Because 
(£,„)" = (£)". m = 1, 2, 3, is assumed in Eq. (4.32). Eqs. (4.33). 
(7.4) and (7.5) require that (i) (V nm )" be in the domain of £(.v), 
and (ii) ()</?< 1. 

Note that ( F mav ) ■' is proportional to A//Ax. Thus. Eqs. (7.3) 
and (7.5) imply that (s)] is an increasing function of At/ Ax, 
i.e., it decreases as Ar decreases if other parameters are held 
constant. Because numerical dissipation decreases as (£)'' de- 
creases, with other factors being equal, the replacement of a 
constant e with (£)■' has an effect in reducing numerical dissipa- 
tion as Ar decreases. This effect will compensate for the ob- 
served opposite effect on numerical dissipation as At decreases 
with e. Ax, and the total marching time is being held constant. 

Moreover, for a fixed a, IV(x_, x + ; a, 0) — > (x + x.)/2 as 
P 0. This fact, coupled with Eq. (4.37), implies that the 
numerical dissipation introduced as a result of replacing 
Wm)j with (rC); will decrease as j3 decreases. Because 
( r^max )/ is proportional to A z/A.r, with other factors being equal, 
the replacement of Ui‘„ n Y‘ by («"„)" defined in Eq. (7.6), has an 
effect in reducing numerical dissipation as A/ decreases. This 
effect will compensate for the observed opposite effect on 
numerical dissipation as A / decreases with a, (i, A.v, and the 
total marching time is held constant. Note that W,(x _v + ; a) 
is a special case of VV(.r , xr, a, ji) with (3 = I. 

Assuming that a = 1 and b = 0.5, the numerical results 
shown in Figs. 12. 13, and 14 are generated with At = 0.004 
(CFL = 0.88), A / = 0.0004 (CFL = 0.088), and At = 0.0001 
(CFL = 0.022), respectively. Note that the results shown in 
Fig. 1 2 are almost identical to those shown in Fig. 7 which were 
generated assuming the same conditions but using a simpler 
marching scheme. However, the results shown in Fig, 13 are 
far less diffusive than their counterparts shown in Fig. 1 1. One 
can conclude from this comparison and the results shown in 
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Fig. 14 that the current modified Euler solver is capable of 
generating accurate numerical solutions even lor the case with 
a very small CFL. 

In the above modified Euler scheme, (e)" and ft are expressed 
as two special functions of (F max );, respectively. They are only 
two among many possible choices. The investigation of other 
choices is a subject to be studied in the future. 

The most general marching scheme presented in Section 4 
is that defined by Eqs. (4.24) and (4.31). It requires several 
matrix multiplications at each mesh points and, therefore, is 
much more costly. Thus, its use is difficult to justify unless a 
substantial gain in accuracy can be made. How this most general 
marching scheme can be applied wisely is left for a future study. 

This completes the numerical study of the Euler solver. We 
conclude this section with a numerical evaluation of the Navier- 
Stokes marching scheme defined by Eqs. (5.20) and (5.29). 
Again the initial conditions defined in Eq. (7.1) are assumed, 
and the numerical solutions are compared with the exact weak 
solution of the Euler equations at t = 0.2. The numerical results 
shown in Figs. 15—17 are generated assuming A t — 0.004, A.v 
= 0.01, y - 1.4, and Pr - 0.72. The value of the Prandtl 
number used here is that tor air at standard conditions. The 
values of the Re L for these figures are 2000, 6000, and 
10,000, respectively. 

From the results shown in these figures, one concludes that, 
for a high-Reynolds-number flow, the shock can be resolved 


FIG. 13. The Euler solution (b — 0.5, a — 1, A/ — 0.0004, CFL — 0.088). 
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FIG. 15. The Navier-Stokes solution (Re, — 2000). 


FIG. 17. The Navier-Stokes solution (Re / = I0(X)0). 
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FIG. 16. The Navier-Stokes solution (Re, = 6(XX)). 




within one mesh interval by the current Navier-Stokes solver. 
Also the contact discontinuity can be resolved within a few 
mesh intervals. Note that these results are obtained without 
using any ad hoc parameters or techniques. Because the Reyn- 
olds number is inversely proportional to the physical viscosity, 
as expected, numerical overshoots and oscillations shown in 
these figures increase slightly as the values of the Reynolds 
number increase. 

Furthermore, through repeated numerical experiments using 
different physical and mesh parameters, it is established that 
the current Navier-Stokes solver is stable if, for all ( /, n ) E (2, 

0<Re,, 0<Pr, (F max )'' < 1 . (7.7) 

However, because a Navier-Stokes problem is fundamentally 
an initial-value/boundary-value problem, the current explicit 
marching scheme obviously cannot model such a problem un- 
less the boundary effect is small, i.e., when the contribution of 
the viscous terms to Eqs. (5.20) and (5.29) is small compared 
to that of the convection terms. In general, this implies that the 
current scheme is applicable only to high-Reynolds-number 
flows. Note that the Leapfrog/Dufort-Frankel and the a-p 
schemes [1] also encounter a similar limitation in modelling 
Eq. (2.1). 

Finally, note that the current Navier-Stokes solver with 
Re, - oc (j.e., the physical viscosity vanishes) and Pr = 0 can 
be considered as a nonlinear extension of the inviseid a-p 
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scheme. Because the latter scheme is neutrally stable, generally 
one would expect that a nonlinear extension of such a scheme 
is unstable. However, it has been shown numerically that the 
current Navier— Stokes solver is stable even tor the above lim- 
iting case as long as (XJ'I < 1 all (y\ w) ^ 

8. CONCLUSIONS AND DISCUSSIONS 

Several key limitations of the finite difference, finite volume, 
finite element, and spectral methods were discussed in Section 
1 . The method of space-time conservation element and solution 
element was conceived to overcome these limitations. 

Using the a- pc scheme as an example, major differences 
between the present method and those mentioned above were 
explained in Section 2. This explicit scheme has the unusual 
property that its stability is limited only by the CFL condition, 
i.e., it is independent of /x. Also, it was shown that its amplifica- 
tion factors are identical to those ot the Leapfrog scheme, it 
jjl = 0, and to those of the DuFort-Frankel scheme, if a = 
0. These coincidences are rather unexpected because the a- 
fx scheme and the above classical schemes are derived from 
completely different perspectives, and the current scheme does 
not reduce to the above classical schemes in the limiting cases. 

The inviscid a-pi scheme is reversible in time. Obviously the 
Euler extension of such a scheme cannot model a physical 
problem that is irreversible in time, e.g., an inviscid flow prob- 
lem involving shocks. Thus, the inviscid version was modified 
in Section 3 to form the a-n scheme. This new scheme has the 
unusual property that numerical dissipation is controlled by an 
adjustable parameter e. As a matter of fact, for all wavelengths, 
numerical dissipation can be uniformly bounded from above 
by an arbitrary small number by choosing a small enough e. 
Stability of the ci-f scheme is limited by the CFL condition 
and 0 < e < 1 . Moreover, if e = 0, the amplification factors 
of the a-e scheme are identical to those of the Leapfrog scheme, 
which has no numerical dissipation. On the other hand, it e = 

1, they unexpectedly become identical to each other and to the 
amplification factor of the highly diffusive Lax scheme. Note 
that, because the Lax scheme is very diffusive and uses a mesh 
that is staggered in time, a two-level scheme using such a mesh 
is often associated with a highly diffusive scheme. The ci-f 
scheme, which also uses a mesh staggered in time, demonstrates 
that such a scheme could be free from numerical dissipation. 

In Section 4, the as scheme was extended to become an 
Euler solver. This solver has the unusual property that numeiicnl 
dissipation at any mesh point (/, n) can be controlled by a set 
of local parameters m = L 2, 3. As in the ci-f scheme, 
stability of the Euler solver is limited by the CFL condition 
and the requirement that, tor all ( /, n), 0 ^ (£„,))■ — l ^ — 

2, 3. Note that an Euler solver using a mesh staggered in time 
is usually highly diffusive for a small CFL number. It was 
shown in Section 7 that the current solver is an exception. It 
can generate highly accurate shock tube solutions with the CFL 
number ranging from 0.88 to 0.022. 


In Section 3, the c/-/x scheme was extended to become a 
Navier-Stokes solver. Stability of this explicit solver is also 
limited only by the CFL condition. Despite the fact that it does 
not use (i) any techniques related to the high-resolution upwind 
methods, and (ii) any ad hoc parameter, it was shown in Section 
7 that the current solver is capable of generating highly accurate 
shock tube solutions. Particularly, shock discontinuites can be 
resolved within one mesh interval. 

A summary of the key results ol the present work has been 
given. Behind these results is a continuous effort to maintain 
the simplicity, generality, and accuracy ot the present method. 
This effort is summarized in the following remarks: 

(a) Simplicity . The current numerical framework rests upon 
only two basic building blocks, i.e., the space-time conservation 
and solution elements. It uses only local discrete variables. 
Also, the set of discrete variables in any one of the numerical 
equations to be solved is associated with a single SE oi a 
few immediately neighboring SEs. Thus, local flexibility is 
preserved and one needs only to deal with a very sparse matrix. 
Moreover, tlux evaluation at an interface separating two CEs 
requires no interpolation or extrapolation. Nor does it require 
the use of an ad hoc flux model. Finally, partly because no 
characleristics-based techniques are used, a numerical scheme 
can be constructed by using only the simplest approximation 
techniques. 

(b) Generality . A guiding principle in the design ot the pres- 
ent method is to limit the use of special assumptions or tech- 
niques that would restrict its use in more general situations. 
Thus we do not use characteristics-based techniques, and we 
try to avoid using ad hoc techniques. 

(c) Accuracy. Because (i) a physical solution ot the conserva- 
tion laws may involve shocks or high-gradient regions, and (ii) 
an accurate numerical simulation ot such a solution is difficult to 
obtain without enforcing flux conservation, the present method 
requires that a numerical solution satisfies (i) the differential 
form of the conservation laws uniformly within an SE, and (ii) 
the integral form over any space-time region that is the union 
of any combination of CEs, In addition, accuracy of the piestnt 
method is aided by treating both (i/j; and (w mi )" as independent 
variables, instead of expressing (»,,„)/ as a finite-difference 
approximation involving (u m Yf s of neighboring mesh points. 
The latter approach may result in poor accuracy in a high- 
gradient region. Also, accuracy is enhanced by the tact that the 
(lux at an interface separating two CEs is evaluated without 
interpolation or extrapolation. Moreover, because flux conser- 
vation is fundamentally a property in space-time, the current 
unified treatment of space and time may also contribute to a 
more accurate simulation of the conservation laws. 

As a result of its simplicity and generality, the current frame- 
work is also very flexible in its ability to generate discretized 
equations such that number of equations can match number of 
unknowns. In |5|, this flexibility is demonstrated in a discussion 
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FIG. 18. A regular space-time mesh. 

on how the current framework can be used to discretize a 2D 
steady incompressible Navier-Stokes problem. In the same 
discussion, the important issue of boundary-condition imple- 
mentation is also addressed. 

Finally, the present paper is concluded with remarks on sev- 
eral extensions of the current basic solvers: 

(a) In [6|, the Euler solver discussed in Sections 4 and 7 was 
extended and applied to more complex flow problems involving 
shock tubes of finite or infinite length. The numerical results 
obtained clearly demonstrate the ability of the extended solver 
to resolve discontinuities accurately even in the presence of 
wave interactions and reflections. 

(b) Several solvers developed in the present paper have been 
extended to solve two-dimensional time-marching problems [7, 
8]. The construction of these extensions are simplified greatly 
by the use of a nontraditional space-time mesh. Its use results 
in the simplest stencil possible, i.e., a tetrahedron in a 3D space- 
time with a vertex at the upper time level and the other three 
at the lower time level. Other discussions of these 2D schemes 
were given in Section 1. 

(c) Extensions to solve 2D steady, incompressible Navier- 
Stokes equations were discussed near the end of Section 2. 


where /, n' = 0, ±1, ±2 The system of equations repre- 

sented by Eq. (A. 1) can be divided into two sets completely 
independent from each other. The first set involves only the 
variables associated with those mesh points marked by dots in 
Fig. 18, and the second set, by crosses. Thus, the solution to 
Eq. (A.l) contains two decoupled solutions. Traditionally the 
von Neumann stability analysis for the Lax scheme is performed 
without taking into account this decoupling nature. Consider a 
solution to Eq. (A.l) in which uf = 1 for all mesh points (/, 
n f ) that are marked by dots, and uf = - 1 for all other (/, n ’). 
In reality, this solution represents the union of two completely 
decoupled constant solutions. However, at any time level, the 
combined solution is represented by a Fourier component of 
the shortest wavelength ( = 2A.v') in the traditional analysis. 
Therefore, two decoupled constant solutions may be wrongly 
perceived as a rapidly-varying solution. For the above reason, 
we shall consider each decoupled solution separately in the 
following von Neumann stability analysis. 

Let n = n'/2,j = f! 2, Ax = 2Ax\ and A t = 2AC. Then 
the mesh depicted in Fig. 18 is identical to that depicted in 
Fig. 2(a) except that those mesh points marked by crosses in 
Fig. 18 have no counterparts in Fig. 2(a). As a result, the 
decoupling nature of Eq. (A.l) will be removed if the Lax 
scheme is expressed using the staggered mesh depicted in Fig. 
2(a), i.e., for all ( j , n) £ II, 


(w)+ 1 Hi + w" 


At/2 


n /2 u\ 
: — + a 


[j t 1/2 


4 j- 1/2 


(A.2 


With the aid of Eq. (2.13), Eq. (A. 2) can be simplified as 

u; = £[(1 + v)u) III + (1 - v)u*+Hh (A.3) 

By applying Eq. (A.3) successively, one has 

= ild + vfu f ; , + 2(1 - 1/2)11; + (1 - I4 2 w ; \,]. (A.4) 

In contrast to Eq. (2.19), Eq. (A.4) implies that w" +l does not 
approach w" as At — » 0. Moreover, by substituting 


Note. To obtain the NASA Technical Memorandums re- 
ferred to in the present paper, please contact the author. 


K" = I G( V. Q) Ye"" (i = V^T, -JT < 8 < n) (A.5) 


APPENDIX A: AN ALTERNATIVE STABILITY ANALYSIS into Eq. (A.4), one concludes that the amplification factor of 
FOR THE LAX AND LEAPFROG/DUFORT-FRANKEL the Lax scheme is given by 
SCHEMES 


With the use of the regular mesh depicted in Fig. 18, the 
Lax scheme for solving Eq. (2.22) can be expressed as 


(»; l , + <_,)/2 , «; l „ 

— r~; f a — — - — ; — = 0, 

At 2 Ax 


(A.l) 


GO, 6) = \cos(0/2) — in sin(0/2)] 2 . (A. 6) 


A comparison among Eqs. (3.12), (3.15), and (A. 6) reveals 
that G { i* = G a - = G(p< 0) when e — L 

Because w' ,rl does not approach u” as At 0. It follows 
from Eq. (A.5) that GO, 0) cannot approach 1 as v — > 0. As 
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a matter of fact. G( v, 0) -» cos’(0/2) as v -► 0. In turn, this 
implies that the Lax scheme is highly diffusive when |t>| is 
small. 

With the use of the regular mesh depicted in Fig. 18, the 
Leapfrog/DuFort-Frankel scheme for solving Eq. (2.1) can be 
expressed as 


K/'m -Hf I 

9 A t' 1 ‘ 2A.v' 

(A. 7) 

+ < , - 1 - _ n 
(57? 

where /, ri = 0, ± 1, ±2, .... Even though Eq. (A.7) is a three- 
level scheme while Eq. (A.l) is a two-level scheme, they have 
the same decoupling nature. The decoupling of Eq. (A.7) can 
be removed if the scheme is expressed with respect to the 
staggered mesh depicted in Fig. 2(a), i.e., for all (j, n) E 12, 


To perform the von Neumann stability analysis for Eq. 
(A. 13), let 

u( 7 \ n) = 0)e ijB (i = V /Z T, -tt < 6 < tt), (A. 14) 

where u*(n, 0) is a 2 X 1 column matrix. Substituting Eq. 
(A. 14) into Eq. (A. 13), one obtains 

u*(/? + 1, 0) = \L(v, £ 0)] 2 u*(h, 6), (A. 15) 

where 

L( v , £ 6) = e m L , + e ,w2 L . (A. 1 6) 

According to Eq. (A. 1 5), [L(v, £ 0 ) | 2 is the amplification matrix. 
Substituting Eq. (A. 12) into Eq. (A. 16), one has 

L(i/,£0) 

2[£cos(0/2) - />sin(0/2)| (1 - €)? 

1 +£ 1 + f 

e m 0 



(A. 17) 


~ u j 

A? 


+ a 


,.n 1/2 _ ,.n 1/2 

1/2 «r 1/2 

Ajc 


^ ' 


1/2 _L ,,/i 1/2 __ ,./i _ ..if 1 

W / +1/2 j 1/2 **; _ 


(A. 8) 


= 0 . 


(A.v/2r 

With the aid of Eq. (2.13), Eq. (A.8) can be simplified as 

( 1 + £)«;' = ( 1 - £)W; ' + (i’+ i)nj:}jl 

( A.V) 

- (v - £)«"* iVf; 

Eq. (A.9) can also be expressed in a two-level form, i.e., 
u(y, n) = L. u(j ~i,n - i) + L u(j + i n - £). (A. 1 0) 

Here 

/ M? \ 

(A. 1 1 ) 


U( '-' ,) \ u n Ul 

\ Uj + 1/2 


for all ( /. // ) E 12 w ith n > 0, and 
v + £ 1 ~ ^ 

L + = | 1 + ? 1 + f |, L = 
0 0 


/- 


-(y- £) 
I +f 


(A. 12) 




By applying Eq. (A. 10) successively, one has 

u(y, n+ 1 ) = (L. ) 2 u(j - 1 , n) + (L.L + L L,)u(j,n) 
+ (L ) : u(./ + 1 . n). 


(A. 13) 


The amplification factors A t given in Eq. (2.21) are the eigen- 
values of the amplification matrix [L(v. £ 0)] 2 . 


APPENDIX B: A SAMPLE PROGRAM FOR SOLVING Sod’s 


C 


C 


5 

c 


c 


SHOCK TUBE PROBLEM 


implicit real *8 (a-h,o-z) 

dimension q(3,1000), qn(3,1000), qx(3 r 1000), qt(3,1000), 

. _ ..... i ..I 


it - 100 
dt - 0.4d-2 
dx - O.ld-l 
ga - 1.4d0 
rhol - l.dO 
ul - o.dO 
pi - l.do 
rhor « 0.125d0 
ur ■ O.dO 
pr *■ o.ldo 
ia = 1 

hdt - dt/ 2 . do 
tt - hdt*df loat ( it) 
qdt - dt/4.d0 
hdx * dx/2.d0 
qdx “ dx/4.dO 
dtx - dt/dx 
al - ga - l.dO 
a2 = 3.d0 - ga 
a3 - a2/2.d0 
a4 * 1 . 5d0*al 
q ( l , 1) = rhol 
q(2,l) - rhol*ul 

q(3,l) - pl/al + 0. 5d0*rhol*ul**2 

itp - it + 1 

do 5 j - 1, itp 

q(l,j+l) - rhor 

q(2,j+l) - rhor*ur 

q(3,j+l) - pr/al + 0 . 5 d 0 *rhor*ur **2 
do 5 i - 1,3 
qx(i, j) ■ O.dO 
continue 

open (unit-8, file-'for008') 
write (8,10) tt r it , ia 
write (8,20) dt,dx,ga 
write (8,30) rhol, ul, pi 
write (8,40) rhor,ur,pr 


m - 2 

do 400 i - 1, it 
do 100 j « l,m 
w2 - q(2, j)/q(l,j) 
W3 - q(3, j)/q(l, j) 
£21 - -a3*w2**2 


£22 “ a2*w2 

£31 - al*w2**3 - ga*w2*v3 
f 32 - ga*w3 - a4*w2**2 
£33 - ga*w2 
qt(l,j) “ -qx( 2 ,j) 

qt ( 2 , j ) - -(£21*qx(l, j) + 
qt ( 3 , j ) - - ( f 3 l*qx ( 1 , j ) 


£2 2 *qx ( 2 , j ) + al*qx(3,j)) 
f 32*qx (2 , j ) + £33*qx(3,j>) 


JC I J . ] ) * -I LJl " 4 * \ J ; T j/ A 

3(1, j) * qdx*qx ( 1 , j ) + dtx*(q(2,j) + qdt*qt ( 2 , j > ) 
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s(2,j) - qdx*qx(2,j) + dtx* ( f 21* (q< 1 , j ) + qdt*qt(l,j)) + 

* f22*(q(2,j) + qdt*qt (2 , j ) ) + al*(q(3,j) + qdt*qt (3 , j) ) ) 
s C 3 , j ) - qdx*qx(3 r j) + dtx* ( f 31* (q ( 1 , j ) + qdt*qt(l,j)) + 

* f32*(q(2,j) + qdt*qt(2,j)> + f33*(q(3,j) + qdt*qt (3 , j ) ) ) 

100 continue 

nan * m - 1 
do 200 j « 1 , bod 
do 200 k =1,3 

qn(k,j+l) = 0. 5d0* (q(k, j) + q(k,j+l) + s(k,j> - s(k,j+l>) 
vxl (k) = (qn (k, j + 1) - q (k, j ) - hdt*qt (k, j > >/hdx 
vxr (k) - (q(k, j+1) + hdt*qt (k, j+1) - qn(k, j+1) )/hdx 
qx{k,j+l) = (vxl {k) * (dabs(vxr(k) ) ) **ia + vxr (k) * (dabs (vxl (k) ) ) 

* **ia)/ ( (dabs (vxl (k) )) **ia + (dabs (vxr (k) )) **ia + l.d-60) 
200 continue 

do 300 j = 2,m 
do 300 k = 1,3 
q(*,j) - qn(k,j) 

300 continue 
m = st + 1 
400 continue 
c 

t2 = dx*df loat (itp) 
xx (1) = -0 . 5d0*t2 
do 500 j = i,itp 
xx (j + 1) = xx(j) + dx 
500 continue 

do 600 j - l,» 
x - q{2, j)/q(l, j) 

z = al* (q ( 3 , j ) - 0 . 5d0*x* *2 *q ( 1 , j ) ) 
write (8,50) xx (j) , q ( l , j ) , x, z 
600 continue 

c 

close (unit=8) 

10 format ( ' t = ',gl4,7,' it = * , i4 , * ia » ',i4) 

20 format ( ' dt = ',gl4.7, # dx = ',gl4.7,' gamma = ',gl4.7) 

30 format ( ' rhol = ',914.7,' ul = ',gl4.7,' pi = ',gl4.7) 

40 format ( ' rhor = ',gl4.7,' ur = ',gl4.7,' pr = ',gl4.7) 

50 format(' x =',f8.4, r rho =',gl4.7,' u =',gl4.7,' p =',gl4.7) 
stop 
end 

400 continue 
c 

t 2 = dx*df loat ( itp) 
xx ( 1 ) = -0 . 5d0*t2 
do 500 j = 1 , itp 
xx (j+ 1 ) = xx(j) + dx 
500 continue 

do 600 j = l,m 
x - q(2, j)/q( 1, j) 

z - al* ( q ( 3 , j ) - 0. 5d0*x**2*q(l, j) ) 
write (8,50) xx { j ) , q (1 , j ) , x, z 
600 continue 
c 

close (unit=8) 

10 format ( ' t = ',gl4.7,' it = ',i4,' ic = # ,i4) 

20 fonnat(' dt = ',gl4.7,' dx = ',gl4.7, # gamma = ',gl4.7) 

30 format ( ' rhol = ',gl4.7,' ul = ',gl4.7,' pi « ',gl4.7) 

40 format ( ' rhor = J ,gl4.7 ( ' ur = ',gl4.7,' pr = ',gl4.7) 

50 format ( ' x =',f8.4,' rho =',gl4.7,' u =',gl4.7,' p =',gl4.7) 

stop 
end 
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